R Workflow Scholier et al.,

Introduction

This Quarto file documents the R worflow used to analyse the data in the paper ‘Effects of past and present habitat on the gut microbiota of a wild rodent’ written by Tiffany Scholier, Anton Lavrinienko, Eva R. Kallio, Phillip C. Watts & Tapio Mappes.

Legend

Throughout the code, the past habitat/site of origin is referred to as sitetr1 and the present habitat/site of transfer is referred to as sitetr2. Urban and rural sites are abbreviated by urb and cont (control), respectively. For example, the term urbcont refers to animals that were translocated from urban forest sites to rural forest sites. Overall, time point 1 refers to pre-transfer samples and time point 3 refers to post-transfer samples.

Loading necessary R packages

Code 
library(phyloseq) #working with phyloseq objects
library(dplyr)    #for %>% function
library(qiime2R)  #for read_qza function
library(ggplot2)  #creating plots
library(ggpubr)   #for ggpaired function
library(vegan)    #for adonis2 function
library(lmerTest) #for lmer function
library(rstatix)  #for add_significance function

packageVersion("phyloseq") #1.44.0
packageVersion("dplyr")    #1.1.3
packageVersion("qiime2R")  #0.99.6
packageVersion("ggplot2")  #3.4.3
packageVersion("ggpubr")   #0.6.0
packageVersion("vegan")    #2.6.4
packageVersion("lmerTest") #3.1.3
packageVersion("rstatix")  #0.7.2

Go to the directory

Code 
setwd("~/Documents/R/RT")

Preparing all necessary R files

Code 
#DATAFRAME FILES

meta<-read.table("metadata_quarto.txt", header=T)
meta$treatment <- factor(meta$treatment, levels = c("contcont", "urbcont", "conturb", "urburb"))
meta.mom<-meta %>% filter(type == "mom")
alpharich<-read_qza("observed_features_vector.qza")$data %>% tibble::rownames_to_column("sampleid")
shannon<-read_qza("shannon_vector.qza")$data %>% tibble::rownames_to_column("sampleid")
faith<-read_qza("faith_pd_vector.qza")$data %>% tibble::rownames_to_column("sampleid")

diff.13<-read.table("differences_1_to_3_from_qiime.txt", header=T)
diff.13$treatment <- factor(diff.13$treatment, levels = c("contcont", "urbcont", "conturb", "urburb"))
diff.13$sitetr1<-as.factor(diff.13$sitetr1)
diff.13$sitetr2<-as.factor(diff.13$sitetr2)
diff.13$change<-as.factor(diff.13$change)

dist.13<-read.table("distances_1_to_3_from_qiime.txt", header=T)
dist.13$treatment <- factor(dist.13$treatment, levels = c("contcont", "conturb", "urbcont", "urburb"))
dist.13$sitetr1<-as.factor(dist.13$sitetr1)
dist.13$sitetr2<-as.factor(dist.13$sitetr2)
dist.13$change<-as.factor(dist.13$change)
dist.13.natives<-subset(dist.13, treatment=="urburb" | treatment == "contcont")

metadata<-meta%>% left_join(alpharich)%>% left_join(shannon)%>% left_join(faith)
metadata$sitetr1<-as.factor(metadata$sitetr1)
metadata$sitetr2<-as.factor(metadata$sitetr2)
metadata$time<-as.factor(metadata$time)
metadata$treatment <- factor(metadata$treatment, levels = c("contcont", "conturb", "urbcont", "urburb"))

metadata.mom<-metadata %>% filter(type == "mom")
metadata.mom.time1 <- metadata.mom %>% filter(time == "1")
metadata.mom.time3 <- metadata.mom %>% filter(time == "3")
metadata.mom.fromurban <- metadata.mom %>% filter(sitetr1 == "urban")
metadata.mom.fromcontrol <- metadata.mom %>% filter(sitetr1 == "control")

metadata.pup<-metadata %>% filter(type == "pup")
metadata.pup$geneticcluster<-as.factor(metadata.pup$geneticcluster)

########################################################################
#RARIFIED PHYLOSEQ

psrt<-qza_to_phyloseq(
  features="rt_table16s_filtermin10freq_rar30671.qza",taxonomy="rt_taxonomysilva.qza",metadata = "metadata_quarto.txt",tree = "rt_rooted-tree.qza")
levels(sample_data(psrt)$sitetr1) <- c("control", "urban")
levels(sample_data(psrt)$sitetr2) <- c("control", "urban")
levels(sample_data(psrt)$treatment) <- c("contcont", "conturban", "urbancont","urbanurban")

psrt.mom = psrt %>% subset_samples(type == "mom")
metadata.psrt.mom <- as(sample_data(psrt.mom), "data.frame")
psrt.mom.time3 = subset_samples(psrt.mom, time != "1")
psrt.mom.time1 = subset_samples(psrt.mom, time == "1")

psrt.pup = psrt %>% subset_samples(type == "pup")
metadata.psrt.pup <- as(sample_data(psrt.pup), "data.frame")

psrt.mom3.pup = subset_samples(psrt, time != "1")
metadata.psrt.mom3.pup <- as(sample_data(psrt.mom3.pup), "data.frame")

psrt.matched = psrt.mom3.pup %>% subset_samples(geneticcluster_2 %in% c("C1", "C5", "C8","C9","C10","C12","C14","C15","C17","C18","C23","C24","C26"))
metadata.matched <- as(sample_data(psrt.matched), "data.frame")

Recreate Figure 1

Experimental design of the reciprocal translocation (RT) experiment on wild bank voles

This figure was made with PowerPoint and GIS software.

Recreate Figure 2

Change in the alpha diversity of the gut microbiota of bank voles in terms of past habitat

Code 
plot.asv.13.urban <- ggpaired(metadata.mom.fromurban, 
                                  x = "time", y = "observed_features", 
                                  order = c("1", "3"),color = "black",fill= "#f00078",
                                  line.color = "gray", line.size = 0.4,id = "ind",
                                  ylab = "ASV richness", xlab = "")+
  ggtitle("Urban")+
  scale_x_discrete(labels=c("1" = "Pre-transfer", "3" = "Post-transfer"))+
  scale_y_continuous(limits = c(150,550))+theme_minimal()

plot.asv.13.rural <- ggpaired(metadata.mom.fromcontrol, 
                              x = "time", y = "observed_features", 
                              order = c("1", "3"),color = "black",fill= "#693aa3",
                              line.color = "gray", line.size = 0.4,id = "ind",
                              ylab = "ASV richness", xlab = "")+
  ggtitle("Rural")+
  scale_x_discrete(labels=c("1" = "Pre-transfer", "3" = "Post-transfer"))+
  scale_y_continuous(limits = c(100,350))+theme_minimal()

plot.asv.t3<-ggplot(data=metadata.mom.time3, aes(x=sitetr1, y=observed_features, fill=sitetr1))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "Post-transfer gut microbiota", x= "",y = "")+
  scale_x_discrete(labels=c("Rural","Urban"))+
  ylab("ASV richness")+
  scale_fill_manual(values=c("#693aa3","#f00078"),name=c(""),
                    labels=c("Rural","Urban"))+
  geom_signif(y_position = c(550), xmin = c(1), xmax = c(2),annotation = c("*"), tip_length = 0)+
  theme_minimal()

plot.change.asv<-ggarrange(plot.asv.13.rural, plot.asv.13.urban,
                           common.legend = FALSE,
                           legend = "right",
                           labels = c("B", "C"),
                           ncol = 1, nrow = 2, align="v"+ theme_void())

plot.asv<-ggarrange(plot.asv.t3, plot.change.asv, 
                           common.legend = FALSE,
                           legend = "none",
                           labels = c("A", ""),
                           ncol = 2, nrow = 1, align="v"+ theme_void())
plot.asv #(6 x 6)

Recreate Figure 3

Post-transfer gut microbiota compositions of bank voles in terms of past and present habitat

Panels were either made with R (A and B) or PowerPoint sotftware (C and D).

Code 
set.seed(123)
bray_dist.mom.time3 = phyloseq::distance(psrt.mom.time3, method="bray")
set.seed(123)
cap.mom.3.bray <- ordinate(physeq = psrt.mom.time3, method = "CAP",distance = bray_dist.mom.time3,formula = ~ sitetr1+sitetr2)
plot.mom.3.bray<-plot_ordination(psrt.mom.time3, cap.mom.3.bray,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Adults")+xlab("CAP Axis.1  [8%]")+ylab("CAP Axis.2  [3.8%]")+theme_minimal()

set.seed(123)
bray_dist.pup.all = phyloseq::distance(psrt.pup, method="bray")
set.seed(123)
cap.pup.all.bray <- ordinate(physeq = psrt.pup, method = "CAP",distance = bray_dist.pup.all,formula = ~ sitetr1+sitetr2)
plot.pup.cap.bray<-plot_ordination(psrt.pup, cap.pup.all.bray,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Offspring")+xlab("CAP Axis.1  [3.9%]")+ylab("CAP Axis.2  [1.9%]")+theme_minimal()

plot.cap.fig3<-ggarrange(plot.mom.3.bray, plot.pup.cap.bray,
                           common.legend = TRUE,
                           legend = "right",
                           labels = c("A", "B"),
                           ncol = 2, nrow = 1, align="v"+ theme_void())
plot.cap.fig3 # 10 x 12

Recreate Table S1

Differences in beta diversity between the pre-transfer gut microbiota of urban and rural bank voles

Code 
#Preparation (calculating distances + generating new metadata file)
set.seed(123)
bray_dist.mom.time1 = phyloseq::distance(psrt.mom.time1, method="bray")
set.seed(123)
jacc_dist.mom.time1 = phyloseq::distance(psrt.mom.time1, method="jaccard")
set.seed(123)
wunifrac_dist.mom.time1 = phyloseq::distance(psrt.mom.time1, method="wunifrac")
set.seed(123)
uunifrac_dist.mom.time1 = phyloseq::distance(psrt.mom.time1, method="uunifrac")
metadata.psrt.time1 <- as(sample_data(psrt.mom.time1), "data.frame")

Bray-Curtis metric

Code 
set.seed(123)
adonis.mom.1<-adonis2(bray_dist.mom.time1~sitetr1, data=metadata.psrt.time1, by="margin",perm=9999)
adonis.mom.1 #R2 = 9%, p-value = 3e-04
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.mom.time1 ~ sitetr1, data = metadata.psrt.time1, permutations = 9999, by = "margin")
         Df SumOfSqs     R2      F Pr(>F)    
sitetr1   1   0.7127 0.0902 2.5776  3e-04 ***
Residual 26   7.1886 0.9098                  
Total    27   7.9013 1.0000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.1 <- betadisper(bray_dist.mom.time1, metadata.psrt.time1$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.1) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value Pr(>F)
Groups     1 0.000690 0.0006905  0.1827 0.6726
Residuals 26 0.098249 0.0037788
Code 
boxplot(mod.mom.1)

Jaccard metric

Code 
set.seed(123)
adonis.mom.1.jacc<-adonis2(jacc_dist.mom.time1~sitetr1, data=metadata.psrt.time1, by="margin", perm=9999)
adonis.mom.1.jacc #R2 = 6.8%, p-value = 3e-04
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.mom.time1 ~ sitetr1, data = metadata.psrt.time1, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)    
sitetr1   1   0.6808 0.06779 1.8907  3e-04 ***
Residual 26   9.3620 0.93221                  
Total    27  10.0428 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.1.jacc <- betadisper(jacc_dist.mom.time1, metadata.psrt.time1$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.1.jacc) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     1 0.000321 0.00032109  0.1865 0.6694
Residuals 26 0.044763 0.00172167
Code 
boxplot(mod.mom.1.jacc)

Weighted Unifrac metric

Code 
set.seed(123)
adonis.mom.1.wu<-adonis2(wunifrac_dist.mom.time1~sitetr1, data=metadata.psrt.time1, by="margin", perm=9999)
adonis.mom.1.wu #p-value > 0.05
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wunifrac_dist.mom.time1 ~ sitetr1, data = metadata.psrt.time1, permutations = 9999, by = "margin")
         Df  SumOfSqs      R2      F Pr(>F)
sitetr1   1 0.0002923 0.05254 1.4418 0.1791
Residual 26 0.0052707 0.94746              
Total    27 0.0055629 1.00000
Code 
set.seed(123)
mod.mom.1.wu <- betadisper(wunifrac_dist.mom.time1, metadata.psrt.time1$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.1.wu) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df     Sum Sq    Mean Sq F value Pr(>F)
Groups     1 0.00001552 1.5517e-05  0.5003 0.4857
Residuals 26 0.00080641 3.1016e-05
Code 
boxplot(mod.mom.1.wu)

Unweighted Unifrac metric

Code 
set.seed(123)
adonis.mom.1.uu<-adonis2(uunifrac_dist.mom.time1~sitetr1, data=metadata.psrt.time1, by="margin", perm=9999)
adonis.mom.1.uu #R2 = 7.7%, p-value = 2e-04
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uunifrac_dist.mom.time1 ~ sitetr1, data = metadata.psrt.time1, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)    
sitetr1   1   0.4871 0.07693 2.1669  2e-04 ***
Residual 26   5.8451 0.92307                  
Total    27   6.3323 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.1.uu <- betadisper(uunifrac_dist.mom.time1, metadata.psrt.time1$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.1.uu) #Heterogeneity!
Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value   Pr(>F)   
Groups     1 0.010445 0.0104451  11.043 0.002651 **
Residuals 26 0.024592 0.0009458                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
boxplot(mod.mom.1.uu) #Higher variance in rural samples in comparison to urban samples

Recreate Table S2

Differences in alpha diversity between the pre-transfer gut microbiota of urban and rural bank voles

Code 
#Test normality
shapiro.test(metadata.mom.time1$observed_features) #p-value = 0.264 = NORMAL DISTRIBUTION
shapiro.test(metadata.mom.time1$shannon_entropy) #p-value = 0.458 = NORMAL DISTRIBUTION
shapiro.test(metadata.mom.time1$faith_pd) #p-value = 0.6019 = NORMAL DISTRIBUTION
#Normality = OK -> use linear models

ASV richness

Code 
set.seed(123)
model_asv_t1<-lm(observed_features~sitetr1, data=metadata.mom.time1)
summary(model_asv_t1)

Call:
lm(formula = observed_features ~ sitetr1, data = metadata.mom.time1)

Residuals:
    Min      1Q  Median      3Q     Max 
-127.83  -61.47  -16.60   52.92  220.62 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   264.833     25.589  10.350 1.03e-10 ***
sitetr1urban    2.542     33.851   0.075    0.941    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 88.64 on 26 degrees of freedom
Multiple R-squared:  0.0002168, Adjusted R-squared:  -0.03824 
F-statistic: 0.005638 on 1 and 26 DF,  p-value: 0.9407

Shannon diversity

Code 
set.seed(123)
model_sh_t1<-lm(shannon_entropy~sitetr1, data=metadata.mom.time1)
summary(model_sh_t1)

Call:
lm(formula = shannon_entropy ~ sitetr1, data = metadata.mom.time1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.45787 -0.46105 -0.03178  0.43075  1.66877 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.3367     0.2283  23.377   <2e-16 ***
sitetr1urban   0.1624     0.3020   0.538    0.595    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7908 on 26 degrees of freedom
Multiple R-squared:  0.011, Adjusted R-squared:  -0.02704 
F-statistic: 0.2891 on 1 and 26 DF,  p-value: 0.5954

Faith’s phylogenetic diversity

Code 
set.seed(123)
model_faith_t1<-lm(faith_pd~sitetr1, data=metadata.mom.time1)
summary(model_faith_t1)

Call:
lm(formula = faith_pd ~ sitetr1, data = metadata.mom.time1)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2540 -2.4222  0.4342  1.6654  5.9643 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   18.2507     0.8940  20.416   <2e-16 ***
sitetr1urban  -0.1107     1.1826  -0.094    0.926    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.097 on 26 degrees of freedom
Multiple R-squared:  0.0003366, Adjusted R-squared:  -0.03811 
F-statistic: 0.008756 on 1 and 26 DF,  p-value: 0.9262

Recreate Table S3

Temporal change in the alpha diversity of the gut microbiota of bank voles within each experimental group

Code 
#Order the dataframes (same individuals should we listed together)
metadata.mom.1.3 <- metadata.mom %>% arrange(ind)
Code 
# Test normality of the differences
shapiro.test(diff.13$asv_1_to_3) #p-value = 0.8119 = NORMAL DISTRIBUTION
shapiro.test(diff.13$shannon_1_to_3) #p-value = 0.9922 = NORMAL DISTRIBUTION
shapiro.test(diff.13$faith_1_to_3) #p-value = 0.5452 = NORMAL DISTRIBUTION

ASV richness

Code 
#Normality = OK -> Use t-tests
stat.test.13<-metadata.mom.1.3 %>%
  group_by(treatment) %>%
  t_test(observed_features ~ time, paired=TRUE,alternative = "two.sided") %>%
  adjust_pvalue(method = "fdr") %>% add_significance("p.adj") 
stat.test.13
# A tibble: 4 × 11
  treatment .y.            group1 group2    n1    n2 statistic    df     p p.adj
  <fct>     <chr>          <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl>
1 contcont  observed_feat… 1      3          7     7     0.865     6 0.42  0.539
2 conturb   observed_feat… 1      3          5     5     0.670     4 0.539 0.539
3 urbcont   observed_feat… 1      3          7     7    -1.01      6 0.349 0.539
4 urburb    observed_feat… 1      3          9     9    -0.749     8 0.476 0.539
# ℹ 1 more variable: p.adj.signif <chr>
Code 
bp <- ggpaired(metadata.mom.1.3, x = "time", y = "observed_features", order = c("1", "3"), color="treatment",
               line.color = "gray", line.size = 0.4,id = "ind",ylab = "ASV richness", xlab = "")+ 
    scale_x_discrete(labels=c("1" = "Pre-transfer", "3" = "Post-transfer"))+
    facet_wrap(~ treatment, nrow=2,ncol=2)
bp

Shannon diversity

Code 
stat.test.13.shannon<-metadata.mom.1.3 %>%
  group_by(treatment) %>%
  t_test(shannon_entropy ~ time, paired=TRUE,alternative = "two.sided") %>%
  adjust_pvalue(method = "fdr") %>%add_significance("p.adj") 
stat.test.13.shannon
# A tibble: 4 × 11
  treatment .y.            group1 group2    n1    n2 statistic    df     p p.adj
  <fct>     <chr>          <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl>
1 contcont  shannon_entro… 1      3          7     7     0.941     6 0.383 0.424
2 conturb   shannon_entro… 1      3          5     5     1.74      4 0.157 0.424
3 urbcont   shannon_entro… 1      3          7     7    -1.26      6 0.253 0.424
4 urburb    shannon_entro… 1      3          9     9     0.842     8 0.424 0.424
# ℹ 1 more variable: p.adj.signif <chr>
Code 
bp.sh <- ggpaired(metadata.mom.1.3, x = "time", y = "shannon_entropy", order = c("1", "3"), color="treatment",
                  line.color = "gray", line.size = 0.4,id = "ind",ylab = "Shannon diversity", xlab = "")+ 
    scale_x_discrete(labels=c("1" = "Pre-transfer", "3" = "Post-transfer"))+ 
    facet_wrap(~ treatment, nrow=2,ncol=2)
bp.sh

Faith’s phylogenetic diversity

Code 
stat.test.13.faith<-metadata.mom.1.3 %>%
  group_by(treatment) %>%
  t_test(faith_pd ~ time, paired=TRUE,alternative = "two.sided") %>%
  adjust_pvalue(method = "fdr") %>%add_significance("p.adj") 
stat.test.13.faith
# A tibble: 4 × 11
  treatment .y.      group1 group2    n1    n2 statistic    df     p p.adj
  <fct>     <chr>    <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl>
1 contcont  faith_pd 1      3          7     7  -0.00543     6 0.996 0.996
2 conturb   faith_pd 1      3          5     5   0.695       4 0.526 0.807
3 urbcont   faith_pd 1      3          7     7  -0.962       6 0.373 0.807
4 urburb    faith_pd 1      3          9     9  -0.538       8 0.605 0.807
# ℹ 1 more variable: p.adj.signif <chr>
Code 
bp.pd <- ggpaired(metadata.mom.1.3, x = "time", y = "faith_pd", order = c("1", "3"), color="treatment",
                  line.color = "gray", line.size = 0.4,id = "ind",ylab = "Faith's PD", xlab = "")+ 
    scale_x_discrete(labels=c("1" = "Pre-transfer", "3" = "Post-transfer"))+ 
    facet_wrap(~ treatment, nrow=2,ncol=2)
bp.pd

Recreate Table S4

Paired-differences in alpha diversity (temporal alpha diversity change) between the gut microbiota of bank voles belonging to different experimental groups

ASV richness

Code 
#Normality = OK (see 'Recreate Table S2')-> Use linear models
model.13.asv <- lm(asv_1_to_3 ~ sitetr1*sitetr2, data = diff.13)
summary(model.13.asv) 

Call:
lm(formula = asv_1_to_3 ~ sitetr1 * sitetr2, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-368.43  -63.29  -11.10   89.34  262.22 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)
(Intercept)                -28.143     57.451  -0.490    0.629
sitetr1urban               106.571     81.248   1.312    0.202
sitetr2urban                -7.257     89.002  -0.082    0.936
sitetr1urban:sitetr2urban  -31.394    117.427  -0.267    0.791

Residual standard error: 152 on 24 degrees of freedom
Multiple R-squared:  0.09611,   Adjusted R-squared:  -0.01687 
F-statistic: 0.8507 on 3 and 24 DF,  p-value: 0.4799
Code 
model.13.asv.2 <- lm(asv_1_to_3 ~ sitetr1+sitetr2, data = diff.13)
summary(model.13.asv.2) 

Call:
lm(formula = asv_1_to_3 ~ sitetr1 + sitetr2, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-360.91  -66.30  -10.27   94.17  267.09 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)    -20.63      49.17  -0.420    0.678
sitetr1urban    91.54      57.56   1.590    0.124
sitetr2urban   -25.29      56.97  -0.444    0.661

Residual standard error: 149.2 on 25 degrees of freedom
Multiple R-squared:  0.09342,   Adjusted R-squared:  0.02089 
F-statistic: 1.288 on 2 and 25 DF,  p-value: 0.2935
Code 
model.13.asv.3 <- lm(asv_1_to_3 ~ sitetr1, data = diff.13)
summary(model.13.asv.3) 

Call:
lm(formula = asv_1_to_3 ~ sitetr1, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-346.69  -62.08   -4.19   86.92  281.31 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)    -31.17      42.39  -0.735    0.469
sitetr1urban    87.85      56.07   1.567    0.129

Residual standard error: 146.8 on 26 degrees of freedom
Multiple R-squared:  0.08627,   Adjusted R-squared:  0.05113 
F-statistic: 2.455 on 1 and 26 DF,  p-value: 0.1292

Shannon diversity

Code 
model.13.shannon <- lm(shannon_1_to_3 ~ sitetr1*sitetr2, data = diff.13)
summary(model.13.shannon)

Call:
lm(formula = shannon_1_to_3 ~ sitetr1 * sitetr2, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8977 -0.6074 -0.1111  0.9135  2.3347 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)
(Intercept)                -0.3849     0.5336  -0.721    0.478
sitetr1urban                1.1707     0.7546   1.551    0.134
sitetr2urban               -0.4004     0.8267  -0.484    0.633
sitetr1urban:sitetr2urban  -0.8348     1.0907  -0.765    0.451

Residual standard error: 1.412 on 24 degrees of freedom
Multiple R-squared:  0.1629,    Adjusted R-squared:  0.05824 
F-statistic: 1.557 on 3 and 24 DF,  p-value: 0.2258
Code 
model.13.shannon.2 <- lm(shannon_1_to_3 ~ sitetr1+sitetr2, data = diff.13)
summary(model.13.shannon.2)

Call:
lm(formula = shannon_1_to_3 ~ sitetr1 + sitetr2, data = diff.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.05313 -0.72255 -0.08885  0.90926  2.23625 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   -0.1851     0.4615  -0.401    0.692
sitetr1urban   0.7710     0.5403   1.427    0.166
sitetr2urban  -0.8799     0.5348  -1.645    0.112

Residual standard error: 1.4 on 25 degrees of freedom
Multiple R-squared:  0.1424,    Adjusted R-squared:  0.07384 
F-statistic: 2.076 on 2 and 25 DF,  p-value: 0.1465

Faith’s phylogenetic diversity

Code 
model.13.faith <- lm(faith_1_to_3 ~ sitetr1*sitetr2, data = diff.13)
summary(model.13.faith)

Call:
lm(formula = faith_1_to_3 ~ sitetr1 * sitetr2, data = diff.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.8544  -2.7661   0.3567   3.9439   8.2502 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)
(Intercept)                0.007656   1.961622   0.004    0.997
sitetr1urban               2.200009   2.774152   0.793    0.436
sitetr2urban              -1.753781   3.038931  -0.577    0.569
sitetr1urban:sitetr2urban  0.475530   4.009479   0.119    0.907

Residual standard error: 5.19 on 24 degrees of freedom
Multiple R-squared:  0.07039,   Adjusted R-squared:  -0.04581 
F-statistic: 0.6057 on 3 and 24 DF,  p-value: 0.6177
Code 
model.13.faith.2 <- lm(faith_1_to_3 ~ sitetr1+sitetr2, data = diff.13)
summary(model.13.faith.2)

Call:
lm(formula = faith_1_to_3 ~ sitetr1 + sitetr2, data = diff.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.7659  -2.7206   0.2429   3.9894   8.1364 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   -0.1062     1.6767  -0.063    0.950
sitetr1urban   2.4277     1.9630   1.237    0.228
sitetr2urban  -1.4806     1.9429  -0.762    0.453

Residual standard error: 5.087 on 25 degrees of freedom
Multiple R-squared:  0.06984,   Adjusted R-squared:  -0.00457 
F-statistic: 0.9386 on 2 and 25 DF,  p-value: 0.4045

All metrics in terms of change

Code 
model.13.asv.change <- lm(asv_1_to_3 ~ change, data = diff.13)
summary(model.13.asv.change)

Call:
lm(formula = asv_1_to_3 ~ change, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-321.00  -77.06   15.00   64.69  307.00 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    31.00      44.23   0.701    0.490
changesame    -20.94      58.52  -0.358    0.723

Residual standard error: 153.2 on 26 degrees of freedom
Multiple R-squared:  0.0049,    Adjusted R-squared:  -0.03337 
F-statistic: 0.128 on 1 and 26 DF,  p-value: 0.7234
Code 
model.13.shannon.change <- lm(shannon_1_to_3 ~ change, data = diff.13)
summary(model.13.shannon.change)

Call:
lm(formula = shannon_1_to_3 ~ change, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9259 -0.7607 -0.2163  1.0557  2.6910 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.1312     0.4201   0.312    0.757
changesame   -0.5524     0.5557  -0.994    0.329

Residual standard error: 1.455 on 26 degrees of freedom
Multiple R-squared:  0.03661,   Adjusted R-squared:  -0.0004424 
F-statistic: 0.9881 on 1 and 26 DF,  p-value: 0.3294
Code 
model.13.faith.change <- lm(faith_1_to_3 ~ change, data = diff.13)
summary(model.13.faith.change)

Call:
lm(formula = faith_1_to_3 ~ change, data = diff.13)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.600  -2.834   0.551   3.047   9.898 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.56025    1.49293   0.375    0.711
changesame  -0.03411    1.97496  -0.017    0.986

Residual standard error: 5.172 on 26 degrees of freedom
Multiple R-squared:  1.147e-05, Adjusted R-squared:  -0.03845 
F-statistic: 0.0002983 on 1 and 26 DF,  p-value: 0.9864

Recreate Table S5

Differences in alpha diversity between the post-transfer gut microbiota of bank voles belonging to different experimental groups

Code 
#Test Normality
shapiro.test(metadata.mom.time3$observed_features) #p-value = 0.0513 = NORMAL DISTRIBUTION
shapiro.test(metadata.mom.time3$shannon_entropy)   #p-value = 0.7244 = NORMAL DISTRIBUTION
shapiro.test(metadata.mom.time3$faith_pd)          #p-value = 0.9103 = NORMAL DISTRIBUTION
#Normality = OK -> use linear models

ASV richness

Code 
model_asv_t3<-lm(observed_features~sitetr1*sitetr2, data=metadata.mom.time3)
summary(model_asv_t3)

Call:
lm(formula = observed_features ~ sitetr1 * sitetr2, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-152.11  -62.11    3.00   53.91  226.89 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                232.143     37.094   6.258 1.81e-06 ***
sitetr1urban               113.714     52.459   2.168   0.0403 *  
sitetr2urban                 3.657     57.466   0.064   0.9498    
sitetr1urban:sitetr2urban  -42.403     75.819  -0.559   0.5812    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 98.14 on 24 degrees of freedom
Multiple R-squared:  0.2114,    Adjusted R-squared:  0.1129 
F-statistic: 2.145 on 3 and 24 DF,  p-value: 0.1209
Code 
model_asv_t3<-lm(observed_features~sitetr1+sitetr2, data=metadata.mom.time3)
summary(model_asv_t3)

Call:
lm(formula = observed_features ~ sitetr1 + sitetr2, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-160.00  -70.00    3.00   58.13  219.00 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    242.29      31.90   7.595 5.98e-08 ***
sitetr1urban    93.41      37.35   2.501   0.0193 *  
sitetr2urban   -20.70      36.97  -0.560   0.5805    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 96.78 on 25 degrees of freedom
Multiple R-squared:  0.2012,    Adjusted R-squared:  0.1373 
F-statistic: 3.148 on 2 and 25 DF,  p-value: 0.06036
Code 
model_asv_t3<-lm(observed_features~sitetr1, data=metadata.mom.time3)
summary(model_asv_t3) #p-value = 0.02

Call:
lm(formula = observed_features ~ sitetr1, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-169.06  -76.31   13.13   58.73  209.94 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    233.67      27.57   8.476 5.88e-09 ***
sitetr1urban    90.40      36.47   2.479     0.02 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 95.5 on 26 degrees of freedom
Multiple R-squared:  0.1911,    Adjusted R-squared:   0.16 
F-statistic: 6.144 on 1 and 26 DF,  p-value: 0.02

Shannon diversity

Code 
model_shannon_t3<-lm(shannon_entropy~sitetr1*sitetr2, data=metadata.mom.time3)
summary(model_shannon_t3)

Call:
lm(formula = shannon_entropy ~ sitetr1 * sitetr2, data = metadata.mom.time3)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.89724 -0.36935  0.09247  0.59053  2.07895 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 4.8399     0.3876  12.487 5.46e-12 ***
sitetr1urban                1.2432     0.5482   2.268   0.0326 *  
sitetr2urban               -0.1320     0.6005  -0.220   0.8279    
sitetr1urban:sitetr2urban  -0.7446     0.7922  -0.940   0.3567    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.026 on 24 degrees of freedom
Multiple R-squared:  0.2296,    Adjusted R-squared:  0.1333 
F-statistic: 2.384 on 3 and 24 DF,  p-value: 0.09432
Code 
model_shannon_t3<-lm(shannon_entropy~sitetr1+sitetr2, data=metadata.mom.time3)
summary(model_shannon_t3)

Call:
lm(formula = shannon_entropy ~ sitetr1 + sitetr2, data = metadata.mom.time3)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.03586 -0.39449 -0.00697  0.64402  1.94033 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.0182     0.3373  14.879 6.29e-14 ***
sitetr1urban   0.8867     0.3948   2.246   0.0338 *  
sitetr2urban  -0.5597     0.3908  -1.432   0.1645    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.023 on 25 degrees of freedom
Multiple R-squared:  0.2012,    Adjusted R-squared:  0.1373 
F-statistic: 3.149 on 2 and 25 DF,  p-value: 0.06031
Code 
model_shannon_t3<-lm(shannon_entropy~sitetr1, data=metadata.mom.time3)
summary(model_shannon_t3) #p-value = 0.0538

Call:
lm(formula = shannon_entropy ~ sitetr1, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2807 -0.4144  0.1836  0.7031  2.0587 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4.7850     0.3013   15.88  6.7e-15 ***
sitetr1urban   0.8051     0.3985    2.02   0.0538 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.044 on 26 degrees of freedom
Multiple R-squared:  0.1357,    Adjusted R-squared:  0.1024 
F-statistic: 4.081 on 1 and 26 DF,  p-value: 0.05378

Faith’s phylogenetic diversity

Code 
model_faith_t3<-lm(faith_pd~sitetr1*sitetr2, data=metadata.mom.time3)
summary(model_faith_t3)

Call:
lm(formula = faith_pd ~ sitetr1 * sitetr2, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2530 -2.2990 -0.0003  2.5007  5.0912 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                17.6274     1.2717  13.861 5.98e-13 ***
sitetr1urban                2.7887     1.7984   1.551    0.134    
sitetr2urban               -0.2396     1.9701  -0.122    0.904    
sitetr1urban:sitetr2urban  -1.1603     2.5993  -0.446    0.659    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.365 on 24 degrees of freedom
Multiple R-squared:  0.1231,    Adjusted R-squared:  0.01353 
F-statistic: 1.123 on 3 and 24 DF,  p-value: 0.3593
Code 
model_faith_t3<-lm(faith_pd~sitetr1+sitetr2, data=metadata.mom.time3)
summary(model_faith_t3)

Call:
lm(formula = faith_pd ~ sitetr1 + sitetr2, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4691 -2.4390 -0.0463  2.8725  5.1923 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   17.9051     1.0912  16.409 6.79e-15 ***
sitetr1urban   2.2332     1.2775   1.748   0.0927 .  
sitetr2urban  -0.9061     1.2644  -0.717   0.4802    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.31 on 25 degrees of freedom
Multiple R-squared:  0.1159,    Adjusted R-squared:  0.04513 
F-statistic: 1.638 on 2 and 25 DF,  p-value: 0.2145
Code 
model_faith_t3<-lm(faith_pd~sitetr1, data=metadata.mom.time3)
summary(model_faith_t3) #p-value = 0.105

Call:
lm(formula = faith_pd ~ sitetr1, data = metadata.mom.time3)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.8655 -2.4988  0.3974  2.3988  5.7020 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   17.5276     0.9466  18.516   <2e-16 ***
sitetr1urban   2.1011     1.2522   1.678    0.105    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.279 on 26 degrees of freedom
Multiple R-squared:  0.0977,    Adjusted R-squared:  0.06299 
F-statistic: 2.815 on 1 and 26 DF,  p-value: 0.1054

Recreate Table S6

Paired-distances in beta diversity (temporal beta diversity change) between the gut microbiota of bank voles belonging to different experimental groups

Code 
#Test normality
shapiro.test(dist.13$bray_1_to_3) #p-value = 0.4519 = NORMAL DISTRIBUTION
shapiro.test(dist.13$jacc_1_to_3) #p-value = 0.4296 = NORMAL DISTRIBUTION
shapiro.test(dist.13$wunifrac_1_to_3) #p-value = 0.01187 = NOT NORMAL
#This needs to be adjusted with a log transformation
dist.13$log_wunifrac_1_to_3 <- log(dist.13$wunifrac_1_to_3)
shapiro.test(dist.13$log_wunifrac_1_to_3) #p-value = 0.3476 = NORMAL DISTRIBUTION
shapiro.test(dist.13$unwunifrac_1_to_3) #0.07002 = NORMAL DISTRIBUTION

Bray-Curtis metric

Code 
lm.aov.13.bray <- lm(bray_1_to_3 ~ sitetr1*sitetr2, data = dist.13)
summary(lm.aov.13.bray)

Call:
lm(formula = bray_1_to_3 ~ sitetr1 * sitetr2, data = dist.13)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.162441 -0.069969  0.003738  0.058870  0.220408 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                0.61869    0.04236  14.605 1.94e-13 ***
sitetr1urban               0.15324    0.05991   2.558   0.0173 *  
sitetr2urban               0.12538    0.06563   1.911   0.0681 .  
sitetr1urban:sitetr2urban -0.17219    0.08658  -1.989   0.0583 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1121 on 24 degrees of freedom
Multiple R-squared:  0.2352,    Adjusted R-squared:  0.1396 
F-statistic:  2.46 on 3 and 24 DF,  p-value: 0.08717
Code 
lm.aov.13.bray.2 <- lm(bray_1_to_3 ~ sitetr1+sitetr2, data = dist.13)
summary(lm.aov.13.bray.2) 

Call:
lm(formula = bray_1_to_3 ~ sitetr1 + sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.20201 -0.06944  0.01582  0.08990  0.17919 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.65991    0.03907  16.891 3.48e-15 ***
sitetr1urban  0.07081    0.04574   1.548    0.134    
sitetr2urban  0.02646    0.04527   0.585    0.564    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1185 on 25 degrees of freedom
Multiple R-squared:  0.1092,    Adjusted R-squared:  0.0379 
F-statistic: 1.532 on 2 and 25 DF,  p-value: 0.2358
Code 
lm.aov.13.bray.3 <- lm(bray_1_to_3 ~ sitetr1, data = dist.13)
summary(lm.aov.13.bray.3) #p-value = 0.107

Call:
lm(formula = bray_1_to_3 ~ sitetr1, data = dist.13)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.213041 -0.059301  0.002861  0.093445  0.168166 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.67093    0.03378  19.864   <2e-16 ***
sitetr1urban  0.07467    0.04468   1.671    0.107    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.117 on 26 degrees of freedom
Multiple R-squared:  0.09699,   Adjusted R-squared:  0.06226 
F-statistic: 2.793 on 1 and 26 DF,  p-value: 0.1067

Jaccard metric

Code 
lm.aov.13.jacc <- lm(jacc_1_to_3 ~ sitetr1*sitetr2, data = dist.13)
summary(lm.aov.13.jacc) 

Call:
lm(formula = jacc_1_to_3 ~ sitetr1 * sitetr2, data = dist.13)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.064476 -0.034023 -0.004301  0.030581  0.107604 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                0.74599    0.01773  42.065   <2e-16 ***
sitetr1urban               0.01821    0.02508   0.726    0.475    
sitetr2urban               0.02802    0.02747   1.020    0.318    
sitetr1urban:sitetr2urban -0.04311    0.03625  -1.189    0.246    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.04692 on 24 degrees of freedom
Multiple R-squared:  0.05732,   Adjusted R-squared:  -0.06051 
F-statistic: 0.4865 on 3 and 24 DF,  p-value: 0.6949
Code 
lm.aov.13.jacc.2 <- lm(jacc_1_to_3 ~ sitetr1+sitetr2, data = dist.13)
summary(lm.aov.13.jacc.2) 

Call:
lm(formula = jacc_1_to_3 ~ sitetr1 + sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.07250 -0.02955 -0.01030  0.03583  0.09958 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.756305   0.015594  48.499   <2e-16 ***
sitetr1urban -0.002425   0.018257  -0.133    0.895    
sitetr2urban  0.003259   0.018070   0.180    0.858    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.04731 on 25 degrees of freedom
Multiple R-squared:  0.001764,  Adjusted R-squared:  -0.07809 
F-statistic: 0.02209 on 2 and 25 DF,  p-value: 0.9782

Weighted Unifrac metric

Code 
lm.aov.13.wunifrac <- lm(log_wunifrac_1_to_3 ~ sitetr1*sitetr2, data = dist.13)
summary(lm.aov.13.wunifrac) 

Call:
lm(formula = log_wunifrac_1_to_3 ~ sitetr1 * sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46018 -0.25177 -0.02192  0.25763  0.66431 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -1.84601    0.12408 -14.878  1.3e-13 ***
sitetr1urban               0.46568    0.17548   2.654   0.0139 *  
sitetr2urban               0.04571    0.19222   0.238   0.8141    
sitetr1urban:sitetr2urban -0.24713    0.25361  -0.974   0.3396    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3283 on 24 degrees of freedom
Multiple R-squared:  0.264, Adjusted R-squared:  0.172 
F-statistic:  2.87 on 3 and 24 DF,  p-value: 0.05747
Code 
lm.aov.13.wunifrac.2 <- lm(log_wunifrac_1_to_3 ~ sitetr1+sitetr2, data = dist.13)
summary(lm.aov.13.wunifrac.2) 

Call:
lm(formula = log_wunifrac_1_to_3 ~ sitetr1 + sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50619 -0.18498 -0.06516  0.29733  0.61830 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.78685    0.10811 -16.529 5.74e-15 ***
sitetr1urban  0.34737    0.12656   2.745   0.0111 *  
sitetr2urban -0.09626    0.12527  -0.768   0.4494    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.328 on 25 degrees of freedom
Multiple R-squared:  0.2349,    Adjusted R-squared:  0.1737 
F-statistic: 3.837 on 2 and 25 DF,  p-value: 0.03521
Code 
lm.aov.13.wunifrac.3 <- lm(log_wunifrac_1_to_3 ~ sitetr1, data = dist.13)
summary(lm.aov.13.wunifrac.3) #p-value (origin) = 0.0125 

Call:
lm(formula = log_wunifrac_1_to_3 ~ sitetr1, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54830 -0.23072 -0.06091  0.25695  0.57619 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.82696    0.09392 -19.452   <2e-16 ***
sitetr1urban  0.33333    0.12425   2.683   0.0125 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3254 on 26 degrees of freedom
Multiple R-squared:  0.2168,    Adjusted R-squared:  0.1867 
F-statistic: 7.197 on 1 and 26 DF,  p-value: 0.01252

Unweighted Unifrac metric

Code 
lm.aov.13.unwunifrac <- lm(unwunifrac_1_to_3 ~ sitetr1*sitetr2, data = dist.13)
summary(lm.aov.13.unwunifrac)

Call:
lm(formula = unwunifrac_1_to_3 ~ sitetr1 * sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.07833 -0.03882 -0.01391  0.04370  0.11902 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                0.416245   0.022107  18.828 7.02e-16 ***
sitetr1urban               0.027492   0.031264   0.879    0.388    
sitetr2urban               0.017493   0.034248   0.511    0.614    
sitetr1urban:sitetr2urban -0.003048   0.045186  -0.067    0.947    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05849 on 24 degrees of freedom
Multiple R-squared:  0.0808,    Adjusted R-squared:  -0.0341 
F-statistic: 0.7032 on 3 and 24 DF,  p-value: 0.5594
Code 
lm.aov.13.unwunifrac.2 <- lm(unwunifrac_1_to_3 ~ sitetr1+sitetr2, data = dist.13)
summary(lm.aov.13.unwunifrac.2) #p-value > 0.05

Call:
lm(formula = unwunifrac_1_to_3 ~ sitetr1 + sitetr2, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.07730 -0.03938 -0.01456  0.04333  0.12004 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.41697    0.01889  22.070   <2e-16 ***
sitetr1urban  0.02603    0.02212   1.177    0.250    
sitetr2urban  0.01574    0.02189   0.719    0.479    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05731 on 25 degrees of freedom
Multiple R-squared:  0.08062,   Adjusted R-squared:  0.007072 
F-statistic: 1.096 on 2 and 25 DF,  p-value: 0.3497

All metrics in terms of change

Code 
#All metrics in terms of change
lm.aov.13.bray.change <- lm(bray_1_to_3 ~ change, data = dist.13)
summary(lm.aov.13.bray.change) #p-value = 0.0757

Call:
lm(formula = bray_1_to_3 ~ change, data = dist.13)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.220667 -0.065626 -0.007383  0.055750  0.212674 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.76032    0.03341   22.75   <2e-16 ***
changesame  -0.08176    0.04420   -1.85   0.0757 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1157 on 26 degrees of freedom
Multiple R-squared:  0.1163,    Adjusted R-squared:  0.08231 
F-statistic: 3.422 on 1 and 26 DF,  p-value: 0.07574
Code 
lm.aov.13.jacc.change <- lm(jacc_1_to_3 ~ change, data = dist.13)
summary(lm.aov.13.jacc.change) #p-value = 0.245

Call:
lm(formula = jacc_1_to_3 ~ change, data = dist.13)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.063108 -0.033264 -0.004496  0.027075  0.108972 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.76829    0.01305   58.86   <2e-16 ***
changesame  -0.02054    0.01727   -1.19    0.245    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.04522 on 26 degrees of freedom
Multiple R-squared:  0.05163,   Adjusted R-squared:  0.01515 
F-statistic: 1.415 on 1 and 26 DF,  p-value: 0.2449
Code 
lm.aov.13.wu.change <- lm(log_wunifrac_1_to_3 ~ change, data = dist.13)
summary(lm.aov.13.wu.change)   #p-value = 0.311

Call:
lm(formula = log_wunifrac_1_to_3 ~ change, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54166 -0.22631 -0.06207  0.20003  0.77992 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.5553     0.1040 -14.952 2.78e-14 ***
changesame   -0.1420     0.1376  -1.032    0.311    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3603 on 26 degrees of freedom
Multiple R-squared:  0.03937,   Adjusted R-squared:  0.002423 
F-statistic: 1.066 on 1 and 26 DF,  p-value: 0.3115
Code 
lm.aov.13.uu.change <- lm(unwunifrac_1_to_3 ~ change, data = dist.13)
summary(lm.aov.13.uu.change)   #p-value = 0.991

Call:
lm(formula = unwunifrac_1_to_3 ~ change, data = dist.13)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08416 -0.05045 -0.01211  0.03699  0.11565 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.4395709  0.0169202  25.979   <2e-16 ***
changesame  0.0002637  0.0223833   0.012    0.991    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05861 on 26 degrees of freedom
Multiple R-squared:  5.34e-06,  Adjusted R-squared:  -0.03846 
F-statistic: 0.0001388 on 1 and 26 DF,  p-value: 0.9907
Code 
plot(bray_1_to_3~treatment, data=dist.13)

Code 
plot(jacc_1_to_3~treatment, data=dist.13)

Code 
plot(log_wunifrac_1_to_3~treatment, data=dist.13)

Code 
plot(unwunifrac_1_to_3~treatment, data=dist.13)

Recreate Table S7

Differences in beta diversity between the post-transfer gut microbiota of bank voles belonging to different experimental groups

Code 
#Preparation (calculating distances + generating new metadata file)
set.seed(123)
bray_dist.mom.time3 = phyloseq::distance(psrt.mom.time3, method="bray")
set.seed(123)
jacc_dist.mom.time3 = phyloseq::distance(psrt.mom.time3, method="jaccard")
set.seed(123)
wunifrac_dist.mom.time3 = phyloseq::distance(psrt.mom.time3, method="wunifrac")
set.seed(123)
uunifrac_dist.mom.time3 = phyloseq::distance(psrt.mom.time3, method="uunifrac")
metadata.psrt.time3 <- as(sample_data(psrt.mom.time3), "data.frame")

Bray-Curtis metric

Code 
set.seed(123)
adonis.mom.3<-adonis2(bray_dist.mom.time3~sitetr1*sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3          #p-value > 0.05
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.mom.time3 ~ sitetr1 * sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
                Df SumOfSqs      R2      F Pr(>F)
sitetr1:sitetr2  1   0.2730 0.03753 1.0668 0.3515
Residual        24   6.1425 0.84443              
Total           27   7.2741 1.00000
Code 
set.seed(123)
adonis.mom.3.2<-adonis2(bray_dist.mom.time3~sitetr1+sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.2 
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.mom.time3 ~ sitetr1 + sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr1   1   0.4324 0.05944 1.6849 0.0286 *
sitetr2   1   0.4702 0.06464 1.8323 0.0155 *
Residual 25   6.4155 0.88197                
Total    27   7.2741 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
adonis.mom.3.3<-adonis2(bray_dist.mom.time3~sitetr1, data=metadata.psrt.time3, by="margin",perm=9999)
adonis.mom.3.3 #p = 0.078
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.mom.time3 ~ sitetr1, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr1   1   0.3884 0.05339 1.4666  0.078 .
Residual 26   6.8857 0.94661                
Total    27   7.2741 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.bray.1 <- betadisper(bray_dist.mom.time3, metadata.psrt.time3$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.bray.1) #Heterogeneity!
Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups     1 0.018118 0.0181184  6.2832 0.01877 *
Residuals 26 0.074974 0.0028836                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
boxplot(mod.mom.bray.1) #Urban > Rural

Code 
set.seed(123)
adonis.mom.3.4<-adonis2(bray_dist.mom.time3~sitetr2, data=metadata.psrt.time3, by="margin",perm=9999)
adonis.mom.3.4 #SIGNIFICANT p = 0.0381 
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.mom.time3 ~ sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr2   1   0.4262 0.05859 1.6182 0.0381 *
Residual 26   6.8479 0.94141                
Total    27   7.2741 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.bray.2 <- betadisper(bray_dist.mom.time3, metadata.psrt.time3$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.bray.2) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value Pr(>F)
Groups     1 0.001958 0.0019577  0.3982 0.5335
Residuals 26 0.127820 0.0049162

Jaccard metric

Code 
set.seed(123)
adonis.mom.3.jacc<-adonis2(jacc_dist.mom.time3~sitetr1*sitetr2, data=metadata.psrt.time3, by="margin",perm=9999)
adonis.mom.3.jacc     #p-value > 0.05
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.mom.time3 ~ sitetr1 * sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
                Df SumOfSqs      R2      F Pr(>F)
sitetr1:sitetr2  1   0.3649 0.03825 1.0664 0.3011
Residual        24   8.2126 0.86087              
Total           27   9.5398 1.00000
Code 
set.seed(123)
adonis.mom.3.2.jacc<-adonis2(jacc_dist.mom.time3~sitetr1+sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.2.jacc 
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.mom.time3 ~ sitetr1 + sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr1   1   0.4739 0.04967 1.3812 0.0379 *
sitetr2   1   0.5237 0.05490 1.5265 0.0132 *
Residual 25   8.5775 0.89913                
Total    27   9.5398 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
adonis.mom.3.3.jacc<-adonis2(jacc_dist.mom.time3~sitetr1, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.3.jacc #p = 0.0948
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.mom.time3 ~ sitetr1, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2     F Pr(>F)  
sitetr1   1   0.4386 0.04598 1.253 0.0948 .
Residual 26   9.1012 0.95402               
Total    27   9.5398 1.00000               
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.jacc.1 <- betadisper(jacc_dist.mom.time3, metadata.psrt.time3$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.jacc.1) #Heterogeneity!
Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups     1 0.008423 0.0084228  6.3903 0.01789 *
Residuals 26 0.034270 0.0013181                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
boxplot(mod.mom.jacc.1) #Urban > Rural

Code 
set.seed(123)
adonis.mom.3.4.jacc<-adonis2(jacc_dist.mom.time3~sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.4.jacc #SIGNIFICANT p = 0.0334
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.mom.time3 ~ sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs     R2     F Pr(>F)  
sitetr2   1   0.4884 0.0512 1.403 0.0334 *
Residual 26   9.0514 0.9488               
Total    27   9.5398 1.0000               
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.jacc.2 <- betadisper(jacc_dist.mom.time3, metadata.psrt.time3$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.jacc.2) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     1 0.000914 0.00091446  0.3879 0.5388
Residuals 26 0.061298 0.00235761

Weighted Unifrac metric

Code 
set.seed(123)
adonis.mom.3.wunifrac<-adonis2(wunifrac_dist.mom.time3~sitetr1*sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.wunifrac #p-value > 0.05
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wunifrac_dist.mom.time3 ~ sitetr1 * sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
                Df  SumOfSqs      R2      F Pr(>F)
sitetr1:sitetr2  1 0.0000751 0.01778 0.4977 0.8134
Residual        24 0.0036202 0.85733              
Total           27 0.0042226 1.00000
Code 
set.seed(123)
adonis.mom.3.2.wunifrac<-adonis2(wunifrac_dist.mom.time3~sitetr1+sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.2.wunifrac 
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wunifrac_dist.mom.time3 ~ sitetr1 + sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df  SumOfSqs      R2      F Pr(>F)  
sitetr1   1 0.0002938 0.06959 1.9880 0.0970 .
sitetr2   1 0.0002535 0.06004 1.7153 0.1342  
Residual 25 0.0036953 0.87511                
Total    27 0.0042226 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
adonis.mom.3.3.wunifrac<-adonis2(wunifrac_dist.mom.time3~sitetr1, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.3.wunifrac #p = 0.122
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wunifrac_dist.mom.time3 ~ sitetr1, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df  SumOfSqs      R2      F Pr(>F)
sitetr1   1 0.0002738 0.06485 1.8029  0.122
Residual 26 0.0039488 0.93515              
Total    27 0.0042226 1.00000
Code 
set.seed(123)
mod.mom.wunifrac.1 <- betadisper(wunifrac_dist.mom.time3, metadata.psrt.time3$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.wunifrac.1) #Heterogeneity!
Analysis of Variance Table

Response: Distances
          Df     Sum Sq    Mean Sq F value   Pr(>F)   
Groups     1 0.00023471 2.3471e-04  9.9812 0.003984 **
Residuals 26 0.00061139 2.3515e-05                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
boxplot(mod.mom.wunifrac.1) #Urban > Rural

Code 
set.seed(123)
adonis.mom.3.4.wunifrac<-adonis2(wunifrac_dist.mom.time3~sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.4.wunifrac #p =  0.1722
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wunifrac_dist.mom.time3 ~ sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df  SumOfSqs     R2      F Pr(>F)
sitetr2   1 0.0002335 0.0553 1.5219 0.1722
Residual 26 0.0039891 0.9447              
Total    27 0.0042226 1.0000
Code 
set.seed(123)
mod.mom.wunifrac.2 <- betadisper(wunifrac_dist.mom.time3, metadata.psrt.time3$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.wunifrac.2) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df     Sum Sq   Mean Sq F value Pr(>F)
Groups     1 0.00000361 3.612e-06  0.1047 0.7488
Residuals 26 0.00089675 3.449e-05

Unweighted Unifrac metric

Code 
set.seed(123)
adonis.mom.3.uunifrac<-adonis2(uunifrac_dist.mom.time3~sitetr1*sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.uunifrac #p-value > 0.05
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uunifrac_dist.mom.time3 ~ sitetr1 * sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
                Df SumOfSqs      R2      F Pr(>F)
sitetr1:sitetr2  1   0.2452 0.03910 1.0803 0.2642
Residual        24   5.4462 0.86869              
Total           27   6.2694 1.00000
Code 
set.seed(123)
adonis.mom.3.2.uunifrac<-adonis2(uunifrac_dist.mom.time3~sitetr1+sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.2.uunifrac 
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uunifrac_dist.mom.time3 ~ sitetr1 + sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr1   1   0.3004 0.04791 1.3194 0.0379 *
sitetr2   1   0.2801 0.04469 1.2306 0.0818 .
Residual 25   5.6914 0.90780                
Total    27   6.2694 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
adonis.mom.3.3.uunifrac<-adonis2(uunifrac_dist.mom.time3~sitetr1, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.3.uunifrac #p = 0.122
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uunifrac_dist.mom.time3 ~ sitetr1, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr1   1   0.2979 0.04752 1.2971 0.0436 *
Residual 26   5.9715 0.95248                
Total    27   6.2694 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.uunifrac.1 <- betadisper(uunifrac_dist.mom.time3, metadata.psrt.time3$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.uunifrac.1) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df    Sum Sq   Mean Sq F value Pr(>F)
Groups     1 0.0000587 5.875e-05  0.0524 0.8207
Residuals 26 0.0291454 1.121e-03
Code 
set.seed(123)
adonis.mom.3.4.uunifrac<-adonis2(uunifrac_dist.mom.time3~sitetr2, data=metadata.psrt.time3, by="margin", perm=9999)
adonis.mom.3.4.uunifrac #p =  0.1722
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uunifrac_dist.mom.time3 ~ sitetr2, data = metadata.psrt.time3, permutations = 9999, by = "margin")
         Df SumOfSqs      R2      F Pr(>F)  
sitetr2   1   0.2777 0.04429 1.2049 0.0996 .
Residual 26   5.9917 0.95571                
Total    27   6.2694 1.00000                
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
mod.mom.uunifrac.2 <- betadisper(uunifrac_dist.mom.time3, metadata.psrt.time3$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
anova(mod.mom.uunifrac.2) #No heterogeneity
Analysis of Variance Table

Response: Distances
          Df   Sum Sq    Mean Sq F value Pr(>F)
Groups     1 0.000084 0.00008399  0.0619 0.8055
Residuals 26 0.035292 0.00135739

Recreate Table S8

Differences in alpha diversity and beta diversity between the post-transfer gut microbiota of offspring belonging to different experimental groups

The offspring data consists of dependent data points since multiple siblings are included into the same dataset. For analyses concerning alpha diversity, we can simply include genetic cluster (i.e., siblings belong to the same group) as a random effect in lmer models. For beta diversity analyses, there are no straightforward methods to include random effects within adonis2 models. To deal with the dependence in our data, we coded a loop (n iterations = 1000) where 1) one sibling per genetic cluster gets randomly selected (33 individuals are left), 2) the adonis2 output gets computed based upon this reduced dataset and 3) the average value for all output generated over 1000 runs gets calculated.

The same output is generated by the loop when 1000 or 10 iterations are used, so the workflow in this quarto document contains the lowest number of iterations for convenience and efficiency.

Code 
shapiro.test(metadata.pup$observed_features)     #p-value = 0.0003368
metadata.pup$observed_features_log <-log(metadata.pup$observed_features)
shapiro.test(metadata.pup$observed_features_log) #p-value = 0.2536
shapiro.test(metadata.pup$shannon_entropy)       #p-value = 0.1673
shapiro.test(metadata.pup$faith_pd)              #p-value = 0.09069

ASV richness

Code 
set.seed(123)
model3<-lmer(observed_features_log~sitetr1*sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: observed_features_log ~ sitetr1 * sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 56

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.65339 -0.65276 -0.08963  0.69150  2.07201 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 0.04456  0.2111  
 Residual                   0.08312  0.2883  
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
                          Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)                 5.5878     0.1095 26.3258  51.012   <2e-16 ***
sitetr1urban               -0.1432     0.1676 31.4393  -0.854    0.399    
sitetr2urban               -0.1057     0.1429 25.8481  -0.739    0.466    
sitetr1urban:sitetr2urban   0.2952     0.2145 30.3422   1.376    0.179    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r sttr2r
sitetr1urbn -0.653              
sitetr2urbn -0.766  0.501       
sttr1rbn:s2  0.511 -0.781 -0.666
Code 
set.seed(123)
model3.1<-lmer(observed_features_log~sitetr1+sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3.1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: observed_features_log ~ sitetr1 + sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 56.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8166 -0.6516 -0.1510  0.6490  2.0116 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 0.04869  0.2207  
 Residual                   0.08236  0.2870  
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   5.51225    0.09635 29.21737  57.208   <2e-16 ***
sitetr1urban  0.03483    0.10691 31.28944   0.326    0.747    
sitetr2urban  0.02495    0.10895 30.92654   0.229    0.820    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r
sitetr1urbn -0.476       
sitetr2urbn -0.664 -0.042

Shannon diversity

Code 
set.seed(123)
model3.sh<-lmer(shannon_entropy~sitetr1*sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3.sh)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: shannon_entropy ~ sitetr1 * sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 187.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.66593 -0.59457 -0.05249  0.64406  1.96930 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 0.3220   0.5675  
 Residual                   0.5685   0.7540  
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
                          Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                5.33561    0.29100 22.81653  18.336 3.78e-15 ***
sitetr1urban               0.05075    0.44469 27.51506   0.114    0.910    
sitetr2urban               0.04243    0.37976 22.40694   0.112    0.912    
sitetr1urban:sitetr2urban  0.08890    0.56925 26.52594   0.156    0.877    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r sttr2r
sitetr1urbn -0.654              
sitetr2urbn -0.766  0.501       
sttr1rbn:s2  0.511 -0.781 -0.667
Code 
set.seed(123)
model3.sh.1<-lmer(shannon_entropy~sitetr1+sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3.sh.1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: shannon_entropy ~ sitetr1 + sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 188.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.69871 -0.60899 -0.06708  0.66770  1.97191 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 0.3034   0.5508  
 Residual                   0.5675   0.7533  
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)    5.3109     0.2459 25.0715  21.597   <2e-16 ***
sitetr1urban   0.1083     0.2731 27.1263   0.396    0.695    
sitetr2urban   0.0821     0.2783 26.7577   0.295    0.770    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r
sitetr1urbn -0.474       
sitetr2urbn -0.665 -0.043

Faith’s phylogenetic diversity

Code 
set.seed(123)
model3.pd<-lmer(faith_pd~sitetr1*sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3.pd)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: faith_pd ~ sitetr1 * sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 357.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.42725 -0.68204 -0.05007  0.65669  1.91753 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 4.783    2.187   
 Residual                   6.472    2.544   
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
                          Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)                 19.237      1.064 27.747  18.083   <2e-16 ***
sitetr1urban                -2.075      1.616 31.847  -1.284   0.2083    
sitetr2urban                -1.429      1.389 27.483  -1.029   0.3125    
sitetr1urban:sitetr2urban    3.518      2.070 31.137   1.699   0.0993 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r sttr2r
sitetr1urbn -0.658              
sitetr2urbn -0.766  0.504       
sttr1rbn:s2  0.514 -0.780 -0.671
Code 
set.seed(123)
model3.pd.1<-lmer(faith_pd~sitetr1+sitetr2 + (1|geneticcluster), data=metadata.pup)
summary(model3.pd.1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: faith_pd ~ sitetr1 + sitetr2 + (1 | geneticcluster)
   Data: metadata.pup

REML criterion at convergence: 363.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5107 -0.6411 -0.0145  0.6748  1.9166 

Random effects:
 Groups         Name        Variance Std.Dev.
 geneticcluster (Intercept) 5.382    2.320   
 Residual                   6.439    2.538   
Number of obs: 72, groups:  geneticcluster, 33

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  18.32127    0.94626 30.38238  19.362   <2e-16 ***
sitetr1urban  0.04267    1.04628 32.17939   0.041    0.968    
sitetr2urban  0.14748    1.06688 31.87568   0.138    0.891    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) sttr1r
sitetr1urbn -0.482       
sitetr2urbn -0.661 -0.041

Bray-Curtis metric

Code 
#################################################
#BRAYCURTIS 1000 iterations = results look the same
#################################################
#INTERACTION = NS (0.4787)
#SOM IN TERMS OF SITETR1 = NS (0.3909), IN TERMS OF SITETR2 = almost significant (0.0612)
#SITETR1 = NS = 0.3926
#SITETR2 = almost significant (0.0577)
#################################################
#adonis2
pvalue.adonis2.interaction<-0
pvalue.adonis2.som.sitetr1<-0
pvalue.adonis2.som.sitetr2<-0
pvalue.bray.sitetr1<-0
pvalue.bray.sitetr2<-0
sumofsqs.1<-0
sumofsqs.2<-0
sumofsqs.res<-0
sumofsqs.tot<-0
sumofsqs.sitetr1.1<-0
sumofsqs.sitetr1.res<-0
sumofsqs.sitetr2.1<-0
sumofsqs.sitetr2.res<-0
r2.1<-0
r2.2<-0
r2.res<-0
r2.sitetr1<-0
r2.sitetr2<-0
f.1<-0
f.2<-0
f.sitetr1<-0
f.sitetr2<-0
#betadisper
pvalue.betadisper.sitetr1 <-0
pvalue.betadisper.sitetr2 <-0

n_iterations <- 10
for (i in 1:n_iterations)  
{
  #Take 1 random sample of pup per mom
  meta.pup.1sibling<-metadata.pup %>% group_by(geneticcluster) %>% dplyr::slice_sample(n = 1) %>% as.data.frame()
  #make this subset into a dataframe
  meta_subset_bray<-as.data.frame(meta.pup.1sibling)
  write.table(meta_subset_bray,"meta_subset_bray.txt",sep="\t",row.names=FALSE, quote = FALSE)
  #make a new phyloseq with only 33 pups instead of 72
  psrt.pup.nosibling<-qza_to_phyloseq(
    features="rt_table16s_filtermin10freq_rar30671.qza",
    taxonomy="rt_taxonomysilva.qza",
    metadata = "meta_subset_bray.txt",
    tree = "rt_rooted-tree.qza"
  )
  metadata.psrt.pup.nosibling <- as(sample_data(psrt.pup.nosibling), "data.frame")
  #Calculate distance for this phyloseq
  set.seed(123)
  bray_dist.pup.nosibling = phyloseq::distance(psrt.pup.nosibling, method="bray")
  
  #Calculate adonis for this phyloseq
  set.seed(123)
  pmod.pup.nosibling_adonis2_1<-adonis2(bray_dist.pup.nosibling~sitetr1+sitetr2+sitetr1:sitetr2, data=metadata.psrt.pup.nosibling, by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling_adonis2_2<-adonis2(bray_dist.pup.nosibling~sitetr1+sitetr2, data=metadata.psrt.pup.nosibling, by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr1<-adonis2(bray_dist.pup.nosibling~sitetr1, data=metadata.psrt.pup.nosibling, by="terms", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr2<-adonis2(bray_dist.pup.nosibling~sitetr2, data=metadata.psrt.pup.nosibling, by="terms", perm=9999)
  
  #calculate group dispersion
  set.seed(123)
  betadisper.pup.nosibling.sitetr1 <- betadisper(bray_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling)$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr1.anova<-anova(betadisper.pup.nosibling.sitetr1)
  set.seed(123)
  betadisper.pup.nosibling.sitetr2 <- betadisper(bray_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling)$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr2.anova<-anova(betadisper.pup.nosibling.sitetr2)
  
  #summate all the pvalues into one object
  
  #INTERACTION: sitetr1 * sitetr2
  pvalue.adonis2.interaction<-pvalue.adonis2.interaction+pmod.pup.nosibling_adonis2_1$`Pr(>F)`[1]
  #SUM: sitetr1 + sitetr2
  pvalue.adonis2.som.sitetr1<-pvalue.adonis2.som.sitetr1+pmod.pup.nosibling_adonis2_2$`Pr(>F)`[1]
  pvalue.adonis2.som.sitetr2<-pvalue.adonis2.som.sitetr2+pmod.pup.nosibling_adonis2_2$`Pr(>F)`[2]
  sumofsqs.1<-sumofsqs.1+pmod.pup.nosibling_adonis2_2$SumOfSqs[1]
  sumofsqs.2<-sumofsqs.2+pmod.pup.nosibling_adonis2_2$SumOfSqs[2]
  sumofsqs.res<-sumofsqs.res+pmod.pup.nosibling_adonis2_2$SumOfSqs[3]
  r2.1<-r2.1+pmod.pup.nosibling_adonis2_2$R2[1]
  r2.2<-r2.2+pmod.pup.nosibling_adonis2_2$R2[2]
  f.1<-f.1+pmod.pup.nosibling_adonis2_2$F[1]
  f.2<-f.2+pmod.pup.nosibling_adonis2_2$F[2]
  #Site of origin: sitetr1
  pvalue.bray.sitetr1<-pvalue.bray.sitetr1 + pmod.pup.nosibling.sitetr1$`Pr(>F)`[1] 
  sumofsqs.sitetr1.1<-sumofsqs.sitetr1.1+pmod.pup.nosibling.sitetr1$SumOfSqs[1]
  sumofsqs.sitetr1.res<-sumofsqs.sitetr1.res+pmod.pup.nosibling.sitetr1$SumOfSqs[2]
  r2.sitetr1<-r2.sitetr1+pmod.pup.nosibling.sitetr1$R2[1]
  f.sitetr1<-f.sitetr1+pmod.pup.nosibling.sitetr1$F[1]
  #Site of transfer: sitetr2
  pvalue.bray.sitetr2<-pvalue.bray.sitetr2 + pmod.pup.nosibling.sitetr2$`Pr(>F)`[1] 
  sumofsqs.sitetr2.1<-sumofsqs.sitetr2.1+pmod.pup.nosibling.sitetr2$SumOfSqs[1]
  sumofsqs.sitetr2.res<-sumofsqs.sitetr2.res+pmod.pup.nosibling.sitetr2$SumOfSqs[2]
  r2.sitetr2<-r2.sitetr2+pmod.pup.nosibling.sitetr2$R2[1]
  f.sitetr2<-f.sitetr2+pmod.pup.nosibling.sitetr2$F[1]
  #Group dispersion
  pvalue.betadisper.sitetr1<-pvalue.betadisper.sitetr1 + betadisper.pup.nosibling.sitetr1.anova$`Pr(>F)`[1]
  pvalue.betadisper.sitetr2<-pvalue.betadisper.sitetr2 + betadisper.pup.nosibling.sitetr2.anova$`Pr(>F)`[1]
  
}
#RESULTS

#INTERACTION: sitetr1*sitetr2
final_pvalue.bray.interaction<-pvalue.adonis2.interaction/n_iterations #p=0.4787
#SOM: sitetr1 + sitetr2
final_pvalue.adonis2.som.sitetr1 <- pvalue.adonis2.som.sitetr1/n_iterations  #0.3909
final_pvalue.adonis2.som.sitetr2 <- pvalue.adonis2.som.sitetr2/n_iterations  #0.0612
final_ss.som.sitetr1<-sumofsqs.1/n_iterations #
final_ss.som.sitetr2<-sumofsqs.2/n_iterations #
final_ss.som.res<-sumofsqs.res/n_iterations 
calculated_ss.som.total <- final_ss.som.sitetr1+final_ss.som.sitetr2+final_ss.som.res
final_r2.som.sitetr1 <- r2.1/n_iterations #0.03211807
final_r2.som.sitetr2 <- r2.2/n_iterations #0.04200681
calculated_r2.som.res = (1-(final_r2.som.sitetr1)-(final_r2.som.sitetr2))
final_f.1 <- f.1/n_iterations 
final_f.2 <- f.2/n_iterations 

#SITE OF ORIGIN
final_pvalue.bray.sitetr1 <- pvalue.bray.sitetr1/n_iterations
final_sumofsqs.sitetr1.1 <- sumofsqs.sitetr1.1/n_iterations
final_sumofsqs.sitetr1.res <- sumofsqs.sitetr1.res/n_iterations
calculated_ss.sitetr1.total <- final_sumofsqs.sitetr1.1+final_sumofsqs.sitetr1.res
final_r2.sitetr1.1 <- r2.sitetr1/n_iterations
calculated_r2.sitetr1.res = (1-final_r2.sitetr1.1)
final_f.sitetr1 <- f.sitetr1/n_iterations

#SITE OF TRANSFER
final_pvalue.bray.sitetr2 <- pvalue.bray.sitetr2/n_iterations
final_sumofsqs.sitetr2.1 <- sumofsqs.sitetr2.1/n_iterations
final_sumofsqs.sitetr2.res <- sumofsqs.sitetr2.res/n_iterations
calculated_ss.sitetr2.total <- final_sumofsqs.sitetr2.1+final_sumofsqs.sitetr2.res
final_r2.sitetr2.1 <- r2.sitetr2/n_iterations
calculated_r2.sitetr2.res = (1-final_r2.sitetr2.1)
final_f.sitetr2 <- f.sitetr2/n_iterations
#GROUP DISPERSION
final_pvalue.betadisper.sitetr1 <- pvalue.betadisper.sitetr1/n_iterations #0.5822573
final_pvalue.betadisper.sitetr2 <- pvalue.betadisper.sitetr2/n_iterations #0.9971147

#MAKE TABLES
df_som_bray <- data.frame(var=c('sitetr1', 'sitetr2', 'Residual', 'Total'),
                          Df=c(1,1,30,32),
                          SumOfSqs=c(final_ss.som.sitetr1,final_ss.som.sitetr2,final_ss.som.res,calculated_ss.som.total),
                          R2=c(final_r2.som.sitetr1,final_r2.som.sitetr2,calculated_r2.som.res,1),
                          "F"=c(final_f.1, final_f.2 ,"",""),
                          "p"=c(final_pvalue.adonis2.som.sitetr1,final_pvalue.adonis2.som.sitetr2,"",""))

df_sitetr1_bray <- data.frame(var=c('sitetr1', 'Residual', 'Total'),
                              Df=c(1,31,32),
                              SumOfSqs=c(final_sumofsqs.sitetr1.1 , final_sumofsqs.sitetr1.res , calculated_ss.sitetr1.total),
                              R2=c(final_r2.som.sitetr1 , calculated_r2.sitetr1.res , "1"),
                              "F"=c(final_f.sitetr1,"",""),
                              "p"=c(final_pvalue.bray.sitetr1,"",""))

df_sitetr2_bray <- data.frame(var=c('sitetr2', 'Residual', 'Total'),
                              Df=c(1,31,32),
                              SumOfSqs=c(final_sumofsqs.sitetr2.1 , final_sumofsqs.sitetr2.res , calculated_ss.sitetr2.total),
                              R2=c(final_r2.som.sitetr2 , calculated_r2.sitetr2.res , "1"),
                              "F"=c(final_f.sitetr2,"",""),
                              "p"=c(final_pvalue.bray.sitetr2,"",""))

Output table treatment

Code 
df_som_bray
       var Df  SumOfSqs         R2                F      p
1  sitetr1  1 0.3077837 0.03211807 1.04107717605331 0.3909
2  sitetr2  1 0.4025464 0.04200681 1.36161134601635 0.0612
3 Residual 30 8.8691907 0.92587512                        
4    Total 32 9.5795208 1.00000000
Code 
write.table(df_som_bray, "adonis2_som_bray.txt")

Output table in terms of site of origin

Code 
df_sitetr1_bray
       var Df  SumOfSqs                 R2                F      p
1  sitetr1  1 0.3111465 0.0321180709382094 1.04031652247225 0.3926
2 Residual 31 9.2717371  0.967531016064008                        
3    Total 32 9.5828836                  1
Code 
write.table(df_sitetr1_bray, "adonis2_sitetr1_bray.txt")

Betadisper p-value site of origin

Code 
final_pvalue.betadisper.sitetr1
[1] 0.5822573

Output table in terms of site of transfer

Code 
df_sitetr2_bray
       var Df  SumOfSqs                 R2                F      p
1  sitetr2  1 0.4059091 0.0420068087242215 1.37116896691662 0.0577
2 Residual 31 9.1769745  0.957642278277996                        
3    Total 32 9.5828836                  1
Code 
write.table(df_sitetr2_bray, "adonis2_sitetr2_bray.txt")

Betadisper p-value site of transfer

Code 
final_pvalue.betadisper.sitetr2
[1] 0.9971147

Jaccard metric

Code 
#################################################
#JACCARD 1000 iterations = results look the same
#################################################
#INTERACTION = NS (p=0.5641)
#SOM IN TERMS OF SITETR1 = NS (0.3782) , IN TERMS OF SITETR2 = almost significant(0.0608)
#SITETR1 = NS (0.3794)
#SITETR2 = almost significant (0.0576)
#################################################
#adonis2
pvalue.adonis2.interaction.jacc<-0
pvalue.adonis2.som.sitetr1.jacc<-0
pvalue.adonis2.som.sitetr2.jacc<-0
pvalue.sitetr1.jacc<-0
pvalue.sitetr2.jacc<-0
sumofsqs.1.jacc<-0
sumofsqs.2.jacc<-0
sumofsqs.res.jacc<-0
sumofsqs.tot.jacc<-0
sumofsqs.sitetr1.1.jacc<-0
sumofsqs.sitetr1.res.jacc<-0
sumofsqs.sitetr2.1.jacc<-0
sumofsqs.sitetr2.res.jacc<-0
r2.1.jacc<-0
r2.2.jacc<-0
r2.res.jacc<-0
r2.sitetr1.jacc<-0
r2.sitetr2.jacc<-0
f.1.jacc<-0
f.2.jacc<-0
f.sitetr1.jacc<-0
f.sitetr2.jacc<-0
#betadisper
pvalue.betadisper.sitetr1.jacc <-0
pvalue.betadisper.sitetr2.jacc <-0

n_iterations <- 10
for (i in 1:n_iterations)  
{
  #Take 1 random sample of pup per mom
  meta.pup.1sibling<-metadata.pup %>% group_by(geneticcluster) %>% dplyr::slice_sample(n = 1) %>% as.data.frame()
  #make this subset into a dataframe
  meta_subset_jacc<-as.data.frame(meta.pup.1sibling)
  write.table(meta_subset_jacc,"meta_subset_jacc.txt",sep="\t",row.names=FALSE, quote = FALSE)
  #make a new phyloseq with only 33 pups instead of 72
  psrt.pup.nosibling.jacc<-qza_to_phyloseq(
    features="rt_table16s_filtermin10freq_rar30671.qza",
    taxonomy="rt_taxonomysilva.qza",
    metadata = "meta_subset_jacc.txt",
    tree = "rt_rooted-tree.qza"
  )
  metadata.psrt.pup.nosibling.jacc <- as(sample_data(psrt.pup.nosibling.jacc), "data.frame")
  #Calculate distance for this phyloseq
  set.seed(123)
  jacc_dist.pup.nosibling = phyloseq::distance(psrt.pup.nosibling.jacc, method="jaccard")
  
  #Calculate adonis for this phyloseq
  set.seed(123)
  pmod.pup.nosibling_adonis2_1.jacc <-adonis2(jacc_dist.pup.nosibling~sitetr1+sitetr2+sitetr1:sitetr2, data=metadata.psrt.pup.nosibling.jacc , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling_adonis2_2.jacc <-adonis2(jacc_dist.pup.nosibling~sitetr1+sitetr2, data=metadata.psrt.pup.nosibling.jacc , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr1.jacc <-adonis2(jacc_dist.pup.nosibling~sitetr1, data=metadata.psrt.pup.nosibling.jacc , by="terms", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr2.jacc <-adonis2(jacc_dist.pup.nosibling~sitetr2, data=metadata.psrt.pup.nosibling.jacc , by="terms", perm=9999)

  #calculate group dispersion
  set.seed(123)
  betadisper.pup.nosibling.sitetr1.jacc  <- betadisper(jacc_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.jacc)$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr1.anova.jacc <-anova(betadisper.pup.nosibling.sitetr1.jacc)
  set.seed(123)
  betadisper.pup.nosibling.sitetr2.jacc  <- betadisper(jacc_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.jacc)$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr2.anova.jacc<-anova(betadisper.pup.nosibling.sitetr2.jacc)
  
  #summate all the pvalues into one object

  #INTERACTION: sitetr1 * sitetr2
  pvalue.adonis2.interaction.jacc<-pvalue.adonis2.interaction.jacc+pmod.pup.nosibling_adonis2_1.jacc$`Pr(>F)`[1]
  #SUM: sitetr1 + sitetr2
  pvalue.adonis2.som.sitetr1.jacc<-pvalue.adonis2.som.sitetr1.jacc+pmod.pup.nosibling_adonis2_2.jacc$`Pr(>F)`[1]
  pvalue.adonis2.som.sitetr2.jacc<-pvalue.adonis2.som.sitetr2.jacc+pmod.pup.nosibling_adonis2_2.jacc$`Pr(>F)`[2]
  sumofsqs.1.jacc<-sumofsqs.1.jacc+pmod.pup.nosibling_adonis2_2.jacc$SumOfSqs[1]
  sumofsqs.2.jacc<-sumofsqs.2.jacc+pmod.pup.nosibling_adonis2_2.jacc$SumOfSqs[2]
  sumofsqs.res.jacc<-sumofsqs.res.jacc+pmod.pup.nosibling_adonis2_2.jacc$SumOfSqs[3]
  r2.1.jacc<-r2.1.jacc+pmod.pup.nosibling_adonis2_2.jacc$R2[1]
  r2.2.jacc<-r2.2.jacc+pmod.pup.nosibling_adonis2_2.jacc$R2[2]
  f.1.jacc<-f.1.jacc+pmod.pup.nosibling_adonis2_2.jacc$F[1]
  f.2.jacc<-f.2.jacc+pmod.pup.nosibling_adonis2_2.jacc$F[2]
  #Site of origin: sitetr1
  pvalue.sitetr1.jacc<-pvalue.sitetr1.jacc + pmod.pup.nosibling.sitetr1.jacc$`Pr(>F)`[1] 
  sumofsqs.sitetr1.1.jacc<-sumofsqs.sitetr1.1.jacc+pmod.pup.nosibling.sitetr1.jacc$SumOfSqs[1]
  sumofsqs.sitetr1.res.jacc<-sumofsqs.sitetr1.res.jacc+pmod.pup.nosibling.sitetr1.jacc$SumOfSqs[2]
  r2.sitetr1.jacc<-r2.sitetr1.jacc+pmod.pup.nosibling.sitetr1.jacc$R2[1]
  f.sitetr1.jacc<-f.sitetr1.jacc+pmod.pup.nosibling.sitetr1.jacc$F[1]
  #Site of transfer: sitetr2
  pvalue.sitetr2.jacc<-pvalue.sitetr2.jacc + pmod.pup.nosibling.sitetr2.jacc$`Pr(>F)`[1] 
  sumofsqs.sitetr2.1.jacc<-sumofsqs.sitetr2.1.jacc+pmod.pup.nosibling.sitetr2.jacc$SumOfSqs[1]
  sumofsqs.sitetr2.res.jacc<-sumofsqs.sitetr2.res.jacc+pmod.pup.nosibling.sitetr2.jacc$SumOfSqs[2]
  r2.sitetr2.jacc<-r2.sitetr2.jacc+pmod.pup.nosibling.sitetr2.jacc$R2[1]
  f.sitetr2.jacc<-f.sitetr2.jacc+pmod.pup.nosibling.sitetr2.jacc$F[1]
  #Group dispersion
  pvalue.betadisper.sitetr1.jacc<-pvalue.betadisper.sitetr1.jacc + betadisper.pup.nosibling.sitetr1.anova.jacc$`Pr(>F)`[1]
  pvalue.betadisper.sitetr2.jacc<-pvalue.betadisper.sitetr2.jacc + betadisper.pup.nosibling.sitetr2.anova.jacc$`Pr(>F)`[1]
  
}
#RESULTS

#INTERACTION: sitetr1*sitetr2
final_pvalue.interaction.jacc<-pvalue.adonis2.interaction.jacc/n_iterations #p=0.5641
#SOM: sitetr1 + sitetr2
final_pvalue.adonis2.som.sitetr1.jacc <- pvalue.adonis2.som.sitetr1.jacc/n_iterations  #0.3782
final_pvalue.adonis2.som.sitetr2.jacc <- pvalue.adonis2.som.sitetr2.jacc/n_iterations  #0.0608
final_ss.som.sitetr1.jacc<-sumofsqs.1.jacc/n_iterations #
final_ss.som.sitetr2.jacc<-sumofsqs.2.jacc/n_iterations #
final_ss.som.res.jacc<-sumofsqs.res.jacc/n_iterations 
calculated_ss.som.total.jacc <- final_ss.som.sitetr1.jacc+final_ss.som.sitetr2.jacc+final_ss.som.res.jacc
final_r2.som.sitetr1.jacc <- r2.1.jacc/n_iterations #0.
final_r2.som.sitetr2.jacc <- r2.2.jacc/n_iterations #0.
calculated_r2.som.res.jacc = (1-(final_r2.som.sitetr1.jacc)-(final_r2.som.sitetr2.jacc))
final_f.1.jacc <- f.1.jacc/n_iterations 
final_f.2.jacc <- f.2.jacc/n_iterations 

#SITE OF ORIGIN
final_pvalue.sitetr1.jacc <- pvalue.sitetr1.jacc/n_iterations #0.3794
final_sumofsqs.sitetr1.1.jacc <- sumofsqs.sitetr1.1.jacc/n_iterations
final_sumofsqs.sitetr1.res.jacc <- sumofsqs.sitetr1.res.jacc/n_iterations
calculated_ss.sitetr1.total.jacc <- final_sumofsqs.sitetr1.1.jacc+final_sumofsqs.sitetr1.res.jacc
final_r2.sitetr1.1.jacc <- r2.sitetr1.jacc/n_iterations
calculated_r2.sitetr1.res.jacc = (1-final_r2.sitetr1.1.jacc)
final_f.sitetr1.jacc <- f.sitetr1.jacc/n_iterations

#SITE OF TRANSFER
final_pvalue.sitetr2.jacc <- pvalue.sitetr2.jacc/n_iterations #0.0576
final_sumofsqs.sitetr2.1.jacc <- sumofsqs.sitetr2.1.jacc/n_iterations
final_sumofsqs.sitetr2.res.jacc <- sumofsqs.sitetr2.res.jacc/n_iterations
calculated_ss.sitetr2.total.jacc <- final_sumofsqs.sitetr2.1.jacc+final_sumofsqs.sitetr2.res.jacc
final_r2.sitetr2.1.jacc <- r2.sitetr2.jacc/n_iterations
calculated_r2.sitetr2.res.jacc = (1-final_r2.sitetr2.1.jacc)
final_f.sitetr2.jacc<- f.sitetr2.jacc/n_iterations

#GROUP DISPERSION
final_pvalue.betadisper.sitetr1.jacc <- pvalue.betadisper.sitetr1.jacc/n_iterations #0.6447998
final_pvalue.betadisper.sitetr2.jacc <- pvalue.betadisper.sitetr2.jacc/n_iterations #0.9854238

#MAKE TABLES
df_som_jaccard <- data.frame(var=c('sitetr1', 'sitetr2', 'Residual', 'Total'),
                          Df=c(1,1,30,32),
                          SumOfSqs=c(final_ss.som.sitetr1.jacc,final_ss.som.sitetr2.jacc,final_ss.som.res.jacc,calculated_ss.som.total.jacc),
                          R2=c(final_r2.som.sitetr1.jacc,final_r2.som.sitetr2.jacc,calculated_r2.som.res.jacc,1),
                          "F"=c(final_f.1, final_f.2.jacc ,"",""),
                          "p"=c(final_pvalue.adonis2.som.sitetr1.jacc,final_pvalue.adonis2.som.sitetr2.jacc,"",""))

df_sitetr1_jaccard <- data.frame(var=c('sitetr1', 'Residual', 'Total'),
                              Df=c(1,31,32),
                              SumOfSqs=c(final_sumofsqs.sitetr1.1.jacc , final_sumofsqs.sitetr1.res.jacc , calculated_ss.sitetr1.total.jacc),
                              R2=c(final_r2.som.sitetr1.jacc , calculated_r2.sitetr1.res.jacc , "1"),
                              "F"=c(final_f.sitetr1.jacc,"",""),
                              "p"=c(final_pvalue.sitetr1.jacc,"",""))

df_sitetr2_jaccard <- data.frame(var=c('sitetr2', 'Residual', 'Total'),
                              Df=c(1,31,32),
                              SumOfSqs=c(final_sumofsqs.sitetr2.1.jacc , final_sumofsqs.sitetr2.res.jacc , calculated_ss.sitetr2.total.jacc),
                              R2=c(final_r2.som.sitetr2.jacc , calculated_r2.sitetr2.res.jacc , "1"),
                              "F"=c(final_f.sitetr2.jacc,"",""),
                              "p"=c(final_pvalue.sitetr2.jacc,"",""))

#################################################

Output table treatment

Code 
df_som_jaccard
       var Df   SumOfSqs         R2                F      p
1  sitetr1  1  0.3856027 0.03194705 1.04107717605331 0.3782
2  sitetr2  1  0.4565100 0.03782169 1.22000684105681 0.0608
3 Residual 30 11.2255916 0.93023126                        
4    Total 32 12.0677043 1.00000000
Code 
write.table(df_som_jaccard, "adonis2_som_jaccard.txt")

Output table in terms of site of origin

Code 
df_sitetr1_jaccard
       var Df   SumOfSqs                R2                F      p
1  sitetr1  1  0.3879553 0.031947051657285 1.02949074806646 0.3794
2 Residual 31 11.6821016  0.96785803570328                        
3    Total 32 12.0700569                 1
Code 
write.table(df_sitetr1_jaccard, "adonis2_sitetr1_jaccard.txt")

Betadisper p-value site of origin

Code 
final_pvalue.betadisper.sitetr1.jacc
[1] 0.6447998

Output table in terms of site of transfer

Code 
df_sitetr2_jaccard
       var Df   SumOfSqs                 R2                F      p
1  sitetr2  1  0.4588626 0.0378216901341921 1.22508838445691 0.0576
2 Residual 31 11.6111944  0.961983397226373                        
3    Total 32 12.0700569                  1
Code 
write.table(df_sitetr2_jaccard, "adonis2_sitetr2_jaccard.txt")

Betadisper p-value site of transfer

Code 
final_pvalue.betadisper.sitetr2.jacc
[1] 0.9854238

Weighted Unifrac metric

Code 
#################################################
#WEIGHTED UNIFRAC 1000 iterations = results look the same
#################################################
#INTERACTION = NS (0.1747)
#SOM IN TERMS OF SITETR1 = NS (0.2498), IN TERMS OF SITETR2 = NS (0.4087)
#SITETR1 = NS (0.2193)
#SITETR2 = NS (0.3679)
#################################################
#adonis2
pvalue.adonis2.interaction.wunifrac<-0
pvalue.adonis2.som.sitetr1.wunifrac<-0
pvalue.adonis2.som.sitetr2.wunifrac<-0
pvalue.sitetr1.wunifrac<-0
pvalue.sitetr2.wunifrac<-0
sumofsqs.1.wunifrac<-0
sumofsqs.2.wunifrac<-0
sumofsqs.res.wunifrac<-0
sumofsqs.tot.wunifrac<-0
sumofsqs.sitetr1.1.wunifrac<-0
sumofsqs.sitetr1.res.wunifrac<-0
sumofsqs.sitetr2.1.wunifrac<-0
sumofsqs.sitetr2.res.wunifrac<-0
r2.1.wunifrac<-0
r2.2.wunifrac<-0
r2.res.wunifrac<-0
r2.sitetr1.wunifrac<-0
r2.sitetr2.wunifrac<-0
f.1.wunifrac<-0
f.2.wunifrac<-0
f.sitetr1.wunifrac<-0
f.sitetr2.wunifrac<-0
#betadisper
pvalue.betadisper.sitetr1.wunifrac <-0
pvalue.betadisper.sitetr2.wunifrac <-0

n_iterations <- 10
for (i in 1:n_iterations)  
{
  #Take 1 random sample of pup per mom
  meta.pup.1sibling<-metadata.pup %>% group_by(geneticcluster) %>% dplyr::slice_sample(n = 1) %>% as.data.frame()
  #make this subset into a dataframe
  meta_subset_wunifrac<-as.data.frame(meta.pup.1sibling)
  write.table(meta_subset_wunifrac,"meta_subset_wunifrac.txt",sep="\t",row.names=FALSE, quote = FALSE)
  #make a new phyloseq with only 33 pups instead of 72
  psrt.pup.nosibling.wunifrac<-qza_to_phyloseq(
    features="rt_table16s_filtermin10freq_rar30671.qza",
    taxonomy="rt_taxonomysilva.qza",
    metadata = "meta_subset_wunifrac.txt",
    tree = "rt_rooted-tree.qza"
  )
  metadata.psrt.pup.nosibling.wunifrac <- as(sample_data(psrt.pup.nosibling.wunifrac), "data.frame")
  #Calculate distance for this phyloseq
  set.seed(123)
  wunifrac_dist.pup.nosibling = phyloseq::distance(psrt.pup.nosibling.wunifrac, method="wunifrac")
  
  #Calculate adonis for this phyloseq
  set.seed(123)
  pmod.pup.nosibling_adonis2_1.wunifrac <-adonis2(wunifrac_dist.pup.nosibling~sitetr1+sitetr2+sitetr1:sitetr2, data=metadata.psrt.pup.nosibling.wunifrac , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling_adonis2_2.wunifrac <-adonis2(wunifrac_dist.pup.nosibling~sitetr1+sitetr2, data=metadata.psrt.pup.nosibling.wunifrac , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr1.wunifrac <-adonis2(wunifrac_dist.pup.nosibling~sitetr1, data=metadata.psrt.pup.nosibling.wunifrac , by="terms", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr2.wunifrac <-adonis2(wunifrac_dist.pup.nosibling~sitetr2, data=metadata.psrt.pup.nosibling.wunifrac , by="terms", perm=9999)
  
  #calculate group dispersion
  set.seed(123)
  betadisper.pup.nosibling.sitetr1.wunifrac  <- betadisper(wunifrac_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.wunifrac)$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr1.anova.wunifrac <-anova(betadisper.pup.nosibling.sitetr1.wunifrac)
  set.seed(123)
  betadisper.pup.nosibling.sitetr2.wunifrac  <- betadisper(wunifrac_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.wunifrac)$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr2.anova.wunifrac<-anova(betadisper.pup.nosibling.sitetr2.wunifrac)
  
  #summate all the pvalues into one object
  
  #INTERACTION: sitetr1 * sitetr2
  pvalue.adonis2.interaction.wunifrac<-pvalue.adonis2.interaction.wunifrac+pmod.pup.nosibling_adonis2_1.wunifrac$`Pr(>F)`[1]
  #SUM: sitetr1 + sitetr2
  pvalue.adonis2.som.sitetr1.wunifrac<-pvalue.adonis2.som.sitetr1.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$`Pr(>F)`[1]
  pvalue.adonis2.som.sitetr2.wunifrac<-pvalue.adonis2.som.sitetr2.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$`Pr(>F)`[2]
  sumofsqs.1.wunifrac<-sumofsqs.1.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$SumOfSqs[1]
  sumofsqs.2.wunifrac<-sumofsqs.2.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$SumOfSqs[2]
  sumofsqs.res.wunifrac<-sumofsqs.res.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$SumOfSqs[3]
  r2.1.wunifrac<-r2.1.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$R2[1]
  r2.2.wunifrac<-r2.2.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$R2[2]
  f.1.wunifrac<-f.1.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$F[1]
  f.2.wunifrac<-f.2.wunifrac+pmod.pup.nosibling_adonis2_2.wunifrac$F[2]
  #Site of origin: sitetr1
  pvalue.sitetr1.wunifrac<-pvalue.sitetr1.wunifrac + pmod.pup.nosibling.sitetr1.wunifrac$`Pr(>F)`[1] 
  sumofsqs.sitetr1.1.wunifrac<-sumofsqs.sitetr1.1.wunifrac+pmod.pup.nosibling.sitetr1.wunifrac$SumOfSqs[1]
  sumofsqs.sitetr1.res.wunifrac<-sumofsqs.sitetr1.res.wunifrac+pmod.pup.nosibling.sitetr1.wunifrac$SumOfSqs[2]
  r2.sitetr1.wunifrac<-r2.sitetr1.wunifrac+pmod.pup.nosibling.sitetr1.wunifrac$R2[1]
  f.sitetr1.wunifrac<-f.sitetr1.wunifrac+pmod.pup.nosibling.sitetr1.wunifrac$F[1]
  #Site of transfer: sitetr2
  pvalue.sitetr2.wunifrac<-pvalue.sitetr2.wunifrac + pmod.pup.nosibling.sitetr2.wunifrac$`Pr(>F)`[1] 
  sumofsqs.sitetr2.1.wunifrac<-sumofsqs.sitetr2.1.wunifrac+pmod.pup.nosibling.sitetr2.wunifrac$SumOfSqs[1]
  sumofsqs.sitetr2.res.wunifrac<-sumofsqs.sitetr2.res.wunifrac+pmod.pup.nosibling.sitetr2.wunifrac$SumOfSqs[2]
  r2.sitetr2.wunifrac<-r2.sitetr2.wunifrac+pmod.pup.nosibling.sitetr2.wunifrac$R2[1]
  f.sitetr2.wunifrac<-f.sitetr2.wunifrac+pmod.pup.nosibling.sitetr2.wunifrac$F[1]
  #Group dispersion
  pvalue.betadisper.sitetr1.wunifrac<-pvalue.betadisper.sitetr1.wunifrac + betadisper.pup.nosibling.sitetr1.anova.wunifrac$`Pr(>F)`[1]
  pvalue.betadisper.sitetr2.wunifrac<-pvalue.betadisper.sitetr2.wunifrac + betadisper.pup.nosibling.sitetr2.anova.wunifrac$`Pr(>F)`[1]
  
}
#RESULTS

#INTERACTION: sitetr1*sitetr2
final_pvalue.interaction.wunifrac<-pvalue.adonis2.interaction.wunifrac/n_iterations #p=0.1747
#SOM: sitetr1 + sitetr2
final_pvalue.adonis2.som.sitetr1.wunifrac <- pvalue.adonis2.som.sitetr1.wunifrac/n_iterations  #0.2498
final_pvalue.adonis2.som.sitetr2.wunifrac <- pvalue.adonis2.som.sitetr2.wunifrac/n_iterations  #0.4087
final_ss.som.sitetr1.wunifrac<-sumofsqs.1.wunifrac/n_iterations #
final_ss.som.sitetr2.wunifrac<-sumofsqs.2.wunifrac/n_iterations #
final_ss.som.res.wunifrac<-sumofsqs.res.wunifrac/n_iterations 
calculated_ss.som.total.wunifrac <- final_ss.som.sitetr1.wunifrac+final_ss.som.sitetr2.wunifrac+final_ss.som.res.wunifrac
final_r2.som.sitetr1.wunifrac <- r2.1.wunifrac/n_iterations #0.
final_r2.som.sitetr2.wunifrac <- r2.2.wunifrac/n_iterations #0.
calculated_r2.som.res.wunifrac = (1-(final_r2.som.sitetr1.wunifrac)-(final_r2.som.sitetr2.wunifrac))
final_f.1.wunifrac <- f.1.wunifrac/n_iterations 
final_f.2.wunifrac <- f.2.wunifrac/n_iterations 

#SITE OF ORIGIN
final_pvalue.sitetr1.wunifrac <- pvalue.sitetr1.wunifrac/n_iterations
final_sumofsqs.sitetr1.1.wunifrac <- sumofsqs.sitetr1.1.wunifrac/n_iterations
final_sumofsqs.sitetr1.res.wunifrac <- sumofsqs.sitetr1.res.wunifrac/n_iterations
calculated_ss.sitetr1.total.wunifrac <- final_sumofsqs.sitetr1.1.wunifrac+final_sumofsqs.sitetr1.res.wunifrac
final_r2.sitetr1.1.wunifrac <- r2.sitetr1.wunifrac/n_iterations
calculated_r2.sitetr1.res.wunifrac = (1-final_r2.sitetr1.1.wunifrac)
final_f.sitetr1.wunifrac <- f.sitetr1.wunifrac/n_iterations

#SITE OF TRANSFER
final_pvalue.sitetr2.wunifrac <- pvalue.sitetr2.wunifrac/n_iterations
final_sumofsqs.sitetr2.1.wunifrac <- sumofsqs.sitetr2.1.wunifrac/n_iterations
final_sumofsqs.sitetr2.res.wunifrac <- sumofsqs.sitetr2.res.wunifrac/n_iterations
calculated_ss.sitetr2.total.wunifrac <- final_sumofsqs.sitetr2.1.wunifrac+final_sumofsqs.sitetr2.res.wunifrac
final_r2.sitetr2.1.wunifrac <- r2.sitetr2.wunifrac/n_iterations
calculated_r2.sitetr2.res.wunifrac = (1-final_r2.sitetr2.1.wunifrac)
final_f.sitetr2.wunifrac<- f.sitetr2.wunifrac/n_iterations

#GROUP DISPERSION
final_pvalue.betadisper.sitetr1.wunifrac <- pvalue.betadisper.sitetr1.wunifrac/n_iterations #0.8
final_pvalue.betadisper.sitetr2.wunifrac <- pvalue.betadisper.sitetr2.wunifrac/n_iterations #0.48

#MAKE TABLES
df_som_wunifrac <- data.frame(var=c('sitetr1', 'sitetr2', 'Residual', 'Total'),
                                 Df=c(1,1,30,32),
                                 SumOfSqs=c(final_ss.som.sitetr1.wunifrac,final_ss.som.sitetr2.wunifrac,final_ss.som.res.wunifrac,calculated_ss.som.total.wunifrac),
                                 R2=c(final_r2.som.sitetr1.wunifrac,final_r2.som.sitetr2.wunifrac,calculated_r2.som.res.wunifrac,1),
                                 "F"=c(final_f.1, final_f.2.wunifrac ,"",""),
                                 "p"=c(final_pvalue.adonis2.som.sitetr1.wunifrac,final_pvalue.adonis2.som.sitetr2.wunifrac,"",""))

df_sitetr1_wunifrac <- data.frame(var=c('sitetr1', 'Residual', 'Total'),
                                     Df=c(1,31,32),
                                     SumOfSqs=c(final_sumofsqs.sitetr1.1.wunifrac , final_sumofsqs.sitetr1.res.wunifrac , calculated_ss.sitetr1.total.wunifrac),
                                     R2=c(final_r2.som.sitetr1.wunifrac , calculated_r2.sitetr1.res.wunifrac , "1"),
                                     "F"=c(final_f.sitetr1.wunifrac,"",""),
                                     "p"=c(final_pvalue.sitetr1.wunifrac,"",""))

df_sitetr2_wunifrac <- data.frame(var=c('sitetr2', 'Residual', 'Total'),
                                     Df=c(1,31,32),
                                     SumOfSqs=c(final_sumofsqs.sitetr2.1.wunifrac , final_sumofsqs.sitetr2.res.wunifrac , calculated_ss.sitetr2.total.wunifrac),
                                     R2=c(final_r2.som.sitetr2.wunifrac , calculated_r2.sitetr2.res.wunifrac , "1"),
                                     "F"=c(final_f.sitetr2.wunifrac,"",""),
                                     "p"=c(final_pvalue.sitetr2.wunifrac,"",""))

#################################################

Output table treatment

Code 
df_som_wunifrac
       var Df     SumOfSqs         R2                F      p
1  sitetr1  1 1.050929e-04 0.03597121 1.04107717605331 0.2498
2  sitetr2  1 9.082028e-05 0.03108598 1.00087076252757 0.4087
3 Residual 30 2.722238e-03 0.93294281                        
4    Total 32 2.918151e-03 1.00000000
Code 
write.table(df_som_wunifrac, "adonis2_som_wunifrac.txt")

Output table in terms of site of origin

Code 
df_sitetr1_wunifrac
       var Df     SumOfSqs                 R2              F      p
1  sitetr1  1 0.0001085249 0.0359712119037808 1.195948661205 0.2193
2 Residual 31 0.0028130582  0.962854063603161                      
3    Total 32 0.0029215832                  1
Code 
write.table(df_sitetr1_wunifrac, "adonis2_sitetr1_wunifrac.txt")

Betadisper p-value site of origin

Code 
final_pvalue.betadisper.sitetr1.wunifrac
[1] 0.8006084

Output table in terms of site of transfer

Code 
df_sitetr2_wunifrac
       var Df     SumOfSqs                R2                F      p
1  sitetr2  1 9.425233e-05 0.031085981043027 1.03342075426691 0.3679
2 Residual 31 2.827331e-03 0.967739294463915                        
3    Total 32 2.921583e-03                 1
Code 
write.table(df_sitetr2_wunifrac, "adonis2_sitetr2_wunifrac.txt")

Betadisper p-value site of transfer

Code 
final_pvalue.betadisper.sitetr2.wunifrac
[1] 0.4862331

Unweighted Unifrac metric

Code 
#################################################
#UNWEIGHTED UNIFRAC 1000 iterations = results look the same
#################################################
#INTERACTION = NS (0.534)
#SOM IN TERMS OF SITETR1 = almost significant (p=0.0845), IN TERMS OF SITETR2 = almost significant (p=0.1037)
#SITETR1 = almost significant (p=0.0876)
#SITETR2 = almost significant (p=0.1073)
#################################################
#################################################
#adonis 2
pvalue.adonis2.interaction.unwunifrac<-0
pvalue.adonis2.som.sitetr1.unwunifrac<-0
pvalue.adonis2.som.sitetr2.unwunifrac<-0
pvalue.sitetr1.unwunifrac<-0
pvalue.sitetr2.unwunifrac<-0
sumofsqs.1.unwunifrac<-0
sumofsqs.2.unwunifrac<-0
sumofsqs.res.unwunifrac<-0
sumofsqs.tot.unwunifrac<-0
sumofsqs.sitetr1.1.unwunifrac<-0
sumofsqs.sitetr1.res.unwunifrac<-0
sumofsqs.sitetr2.1.unwunifrac<-0
sumofsqs.sitetr2.res.unwunifrac<-0
r2.1.unwunifrac<-0
r2.2.unwunifrac<-0
r2.res.unwunifrac<-0
r2.sitetr1.unwunifrac<-0
r2.sitetr2.unwunifrac<-0
f.1.unwunifrac<-0
f.2.unwunifrac<-0
f.sitetr1.unwunifrac<-0
f.sitetr2.unwunifrac<-0
#betadisper
pvalue.betadisper.sitetr1.unwunifrac <-0
pvalue.betadisper.sitetr2.unwunifrac <-0

n_iterations <- 10
for (i in 1:n_iterations)  
{
  #Take 1 random sample of pup per mom
  meta.pup.1sibling<-metadata.pup %>% group_by(geneticcluster) %>% dplyr::slice_sample(n = 1) %>% as.data.frame()
  #make this subset into a dataframe
  meta_subset_unwunifrac<-as.data.frame(meta.pup.1sibling)
  write.table(meta_subset_unwunifrac,"meta_subset_unwunifrac.txt",sep="\t",row.names=FALSE, quote = FALSE)
  #make a new phyloseq with only 33 pups instead of 72
  psrt.pup.nosibling.unwunifrac<-qza_to_phyloseq(
    features="rt_table16s_filtermin10freq_rar30671.qza",
    taxonomy="rt_taxonomysilva.qza",
    metadata = "meta_subset_unwunifrac.txt",
    tree = "rt_rooted-tree.qza"
  )
  metadata.psrt.pup.nosibling.unwunifrac <- as(sample_data(psrt.pup.nosibling.unwunifrac), "data.frame")
  #Calculate distance for this phyloseq
  set.seed(123)
  unwunifrac_dist.pup.nosibling = phyloseq::distance(psrt.pup.nosibling.unwunifrac, method="unwunifrac")
  
  #Calculate adonis for this phyloseq
  set.seed(123)
  pmod.pup.nosibling_adonis2_1.unwunifrac <-adonis2(unwunifrac_dist.pup.nosibling~sitetr1+sitetr2+sitetr1:sitetr2, data=metadata.psrt.pup.nosibling.unwunifrac , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling_adonis2_2.unwunifrac <-adonis2(unwunifrac_dist.pup.nosibling~sitetr1+sitetr2, data=metadata.psrt.pup.nosibling.unwunifrac , by="margin", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr1.unwunifrac <-adonis2(unwunifrac_dist.pup.nosibling~sitetr1, data=metadata.psrt.pup.nosibling.unwunifrac , by="terms", perm=9999)
  set.seed(123)
  pmod.pup.nosibling.sitetr2.unwunifrac <-adonis2(unwunifrac_dist.pup.nosibling~sitetr2, data=metadata.psrt.pup.nosibling.unwunifrac , by="terms", perm=9999)
  
  #calculate group dispersion
  set.seed(123)
  betadisper.pup.nosibling.sitetr1.unwunifrac  <- betadisper(unwunifrac_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.unwunifrac)$sitetr1, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr1.anova.unwunifrac <-anova(betadisper.pup.nosibling.sitetr1.unwunifrac)
  set.seed(123)
  betadisper.pup.nosibling.sitetr2.unwunifrac  <- betadisper(unwunifrac_dist.pup.nosibling, phyloseq::sample_data(psrt.pup.nosibling.unwunifrac)$sitetr2, bias.adjust = TRUE, type=c("centroid"), sqrt.dist = FALSE, add = FALSE)
  betadisper.pup.nosibling.sitetr2.anova.unwunifrac<-anova(betadisper.pup.nosibling.sitetr2.unwunifrac)
  
  #summate all the pvalues into one object
  
  #INTERACTION: sitetr1 * sitetr2
  pvalue.adonis2.interaction.unwunifrac<-pvalue.adonis2.interaction.unwunifrac+pmod.pup.nosibling_adonis2_1.unwunifrac$`Pr(>F)`[1]
  #SUM: sitetr1 + sitetr2
  pvalue.adonis2.som.sitetr1.unwunifrac<-pvalue.adonis2.som.sitetr1.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$`Pr(>F)`[1]
  pvalue.adonis2.som.sitetr2.unwunifrac<-pvalue.adonis2.som.sitetr2.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$`Pr(>F)`[2]
  sumofsqs.1.unwunifrac<-sumofsqs.1.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$SumOfSqs[1]
  sumofsqs.2.unwunifrac<-sumofsqs.2.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$SumOfSqs[2]
  sumofsqs.res.unwunifrac<-sumofsqs.res.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$SumOfSqs[3]
  r2.1.unwunifrac<-r2.1.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$R2[1]
  r2.2.unwunifrac<-r2.2.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$R2[2]
  f.1.unwunifrac<-f.1.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$F[1]
  f.2.unwunifrac<-f.2.unwunifrac+pmod.pup.nosibling_adonis2_2.unwunifrac$F[2]
  #Site of origin: sitetr1
  pvalue.sitetr1.unwunifrac<-pvalue.sitetr1.unwunifrac + pmod.pup.nosibling.sitetr1.unwunifrac$`Pr(>F)`[1] 
  sumofsqs.sitetr1.1.unwunifrac<-sumofsqs.sitetr1.1.unwunifrac+pmod.pup.nosibling.sitetr1.unwunifrac$SumOfSqs[1]
  sumofsqs.sitetr1.res.unwunifrac<-sumofsqs.sitetr1.res.unwunifrac+pmod.pup.nosibling.sitetr1.unwunifrac$SumOfSqs[2]
  r2.sitetr1.unwunifrac<-r2.sitetr1.unwunifrac+pmod.pup.nosibling.sitetr1.unwunifrac$R2[1]
  f.sitetr1.unwunifrac<-f.sitetr1.unwunifrac+pmod.pup.nosibling.sitetr1.unwunifrac$F[1]
  #Site of transfer: sitetr2
  pvalue.sitetr2.unwunifrac<-pvalue.sitetr2.unwunifrac + pmod.pup.nosibling.sitetr2.unwunifrac$`Pr(>F)`[1] 
  sumofsqs.sitetr2.1.unwunifrac<-sumofsqs.sitetr2.1.unwunifrac+pmod.pup.nosibling.sitetr2.unwunifrac$SumOfSqs[1]
  sumofsqs.sitetr2.res.unwunifrac<-sumofsqs.sitetr2.res.unwunifrac+pmod.pup.nosibling.sitetr2.unwunifrac$SumOfSqs[2]
  r2.sitetr2.unwunifrac<-r2.sitetr2.unwunifrac+pmod.pup.nosibling.sitetr2.unwunifrac$R2[1]
  f.sitetr2.unwunifrac<-f.sitetr2.unwunifrac+pmod.pup.nosibling.sitetr2.unwunifrac$F[1]
  #Group dispersion
  pvalue.betadisper.sitetr1.unwunifrac<-pvalue.betadisper.sitetr1.unwunifrac + betadisper.pup.nosibling.sitetr1.anova.unwunifrac$`Pr(>F)`[1]
  pvalue.betadisper.sitetr2.unwunifrac<-pvalue.betadisper.sitetr2.unwunifrac + betadisper.pup.nosibling.sitetr2.anova.unwunifrac$`Pr(>F)`[1]
  
}
#RESULTS

#INTERACTION: sitetr1*sitetr2
final_pvalue.interaction.unwunifrac<-pvalue.adonis2.interaction.unwunifrac/n_iterations #p=0.534
#SOM: sitetr1 + sitetr2
final_pvalue.adonis2.som.sitetr1.unwunifrac <- pvalue.adonis2.som.sitetr1.unwunifrac/n_iterations  #0.0845
final_pvalue.adonis2.som.sitetr2.unwunifrac <- pvalue.adonis2.som.sitetr2.unwunifrac/n_iterations  #0.1037
final_ss.som.sitetr1.unwunifrac<-sumofsqs.1.unwunifrac/n_iterations #
final_ss.som.sitetr2.unwunifrac<-sumofsqs.2.unwunifrac/n_iterations #
final_ss.som.res.unwunifrac<-sumofsqs.res.unwunifrac/n_iterations 
calculated_ss.som.total.unwunifrac <- final_ss.som.sitetr1.unwunifrac+final_ss.som.sitetr2.unwunifrac+final_ss.som.res.unwunifrac
final_r2.som.sitetr1.unwunifrac <- r2.1.unwunifrac/n_iterations #0.
final_r2.som.sitetr2.unwunifrac <- r2.2.unwunifrac/n_iterations #0.
calculated_r2.som.res.unwunifrac = (1-(final_r2.som.sitetr1.unwunifrac)-(final_r2.som.sitetr2.unwunifrac))
final_f.1.unwunifrac <- f.1.unwunifrac/n_iterations 
final_f.2.unwunifrac <- f.2.unwunifrac/n_iterations 

#SITE OF ORIGIN
final_pvalue.sitetr1.unwunifrac <- pvalue.sitetr1.unwunifrac/n_iterations
final_sumofsqs.sitetr1.1.unwunifrac <- sumofsqs.sitetr1.1.unwunifrac/n_iterations
final_sumofsqs.sitetr1.res.unwunifrac <- sumofsqs.sitetr1.res.unwunifrac/n_iterations
calculated_ss.sitetr1.total.unwunifrac <- final_sumofsqs.sitetr1.1.unwunifrac+final_sumofsqs.sitetr1.res.unwunifrac
final_r2.sitetr1.1.unwunifrac <- r2.sitetr1.unwunifrac/n_iterations
calculated_r2.sitetr1.res.unwunifrac = (1-final_r2.sitetr1.1.unwunifrac)
final_f.sitetr1.unwunifrac <- f.sitetr1.unwunifrac/n_iterations

#SITE OF TRANSFER
final_pvalue.sitetr2.unwunifrac <- pvalue.sitetr2.unwunifrac/n_iterations
final_sumofsqs.sitetr2.1.unwunifrac <- sumofsqs.sitetr2.1.unwunifrac/n_iterations
final_sumofsqs.sitetr2.res.unwunifrac <- sumofsqs.sitetr2.res.unwunifrac/n_iterations
calculated_ss.sitetr2.total.unwunifrac <- final_sumofsqs.sitetr2.1.unwunifrac+final_sumofsqs.sitetr2.res.unwunifrac
final_r2.sitetr2.1.unwunifrac <- r2.sitetr2.unwunifrac/n_iterations
calculated_r2.sitetr2.res.unwunifrac = (1-final_r2.sitetr2.1.unwunifrac)
final_f.sitetr2.unwunifrac<- f.sitetr2.unwunifrac/n_iterations

#GROUP DISPERSION
final_pvalue.betadisper.sitetr1.unwunifrac <- pvalue.betadisper.sitetr1.unwunifrac/n_iterations #0.2114821
final_pvalue.betadisper.sitetr2.unwunifrac <- pvalue.betadisper.sitetr2.unwunifrac/n_iterations #0.827527

#MAKE TABLES
df_som_unwunifrac <- data.frame(var=c('sitetr1', 'sitetr2', 'Residual', 'Total'),
                                Df=c(1,1,30,32),
                                SumOfSqs=c(final_ss.som.sitetr1.unwunifrac,final_ss.som.sitetr2.unwunifrac,final_ss.som.res.unwunifrac,calculated_ss.som.total.unwunifrac),
                                R2=c(final_r2.som.sitetr1.unwunifrac,final_r2.som.sitetr2.unwunifrac,calculated_r2.som.res.unwunifrac,1),
                                "F"=c(final_f.1, final_f.2.unwunifrac ,"",""),
                                "p"=c(final_pvalue.adonis2.som.sitetr1.unwunifrac,final_pvalue.adonis2.som.sitetr2.unwunifrac,"",""))

df_sitetr1_unwunifrac <- data.frame(var=c('sitetr1', 'Residual', 'Total'),
                                    Df=c(1,31,32),
                                    SumOfSqs=c(final_sumofsqs.sitetr1.1.unwunifrac , final_sumofsqs.sitetr1.res.unwunifrac , calculated_ss.sitetr1.total.unwunifrac),
                                    R2=c(final_r2.som.sitetr1.unwunifrac , calculated_r2.sitetr1.res.unwunifrac , "1"),
                                    "F"=c(final_f.sitetr1.unwunifrac,"",""),
                                    "p"=c(final_pvalue.sitetr1.unwunifrac,"",""))

df_sitetr2_unwunifrac <- data.frame(var=c('sitetr2', 'Residual', 'Total'),
                                    Df=c(1,31,32),
                                    SumOfSqs=c(final_sumofsqs.sitetr2.1.unwunifrac , final_sumofsqs.sitetr2.res.unwunifrac , calculated_ss.sitetr2.total.unwunifrac),
                                    R2=c(final_r2.som.sitetr2.unwunifrac , calculated_r2.sitetr2.res.unwunifrac , "1"),
                                    "F"=c(final_f.sitetr2.unwunifrac,"",""),
                                    "p"=c(final_pvalue.sitetr2.unwunifrac,"",""))

#################################################

Output table treatment

Code 
df_som_unwunifrac
       var Df  SumOfSqs         R2                F      p
1  sitetr1  1 0.2884181 0.03746472 1.04107717605331 0.0845
2  sitetr2  1 0.2835715 0.03683516 1.19392134628199 0.1037
3 Residual 30 7.1253819 0.92570012                        
4    Total 32 7.6973715 1.00000000
Code 
write.table(df_som_unwunifrac, "adonis2_som_unwunifrac.txt")

Output table in terms of site of origin

Code 
df_sitetr1_unwunifrac
       var Df  SumOfSqs                 R2                F      p
1  sitetr1  1 0.2894394 0.0374647167688244 1.21105103662537 0.0876
2 Residual 31 7.4089534  0.962402622775384                        
3    Total 32 7.6983928                  1
Code 
write.table(df_sitetr1_unwunifrac, "adonis2_sitetr1_unwunifrac.txt")

Betadisper p-value site of origin

Code 
final_pvalue.betadisper.sitetr1.unwunifrac
[1] 0.2114821

Output table in terms of site of transfer

Code 
df_sitetr2_unwunifrac
       var Df  SumOfSqs                 R2                F      p
1  sitetr2  1 0.2845928 0.0368351584366053 1.18999386328529 0.1073
2 Residual 31 7.4138000  0.963032181107603                        
3    Total 32 7.6983928                  1
Code 
write.table(df_sitetr2_unwunifrac, "adonis2_sitetr2_unwunifrac.txt")

Betadisper p-value site of transfer

Code 
final_pvalue.betadisper.sitetr2.unwunifrac
[1] 0.827527

All metrics in terms of genetic cluster

Code 
#ALPHA DIVERSITY
set.seed(123)
model3.geneticcluster<-lm(observed_features_log~geneticcluster, data=metadata.pup)
summary(model3.geneticcluster) #Adjusted R2 = 3256, p-value = 0.01551

Call:
lm(formula = observed_features_log ~ geneticcluster, data = metadata.pup)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.5381 -0.1816  0.0000  0.1154  0.5605 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       5.659482   0.294061  19.246   <2e-16 ***
geneticcluster5  -0.639593   0.339553  -1.884   0.0671 .  
geneticcluster8  -0.342192   0.339553  -1.008   0.3198    
geneticcluster9  -0.408870   0.360150  -1.135   0.2632    
geneticcluster10 -0.221403   0.415865  -0.532   0.5975    
geneticcluster12  0.369037   0.360150   1.025   0.3118    
geneticcluster14 -0.472128   0.360150  -1.311   0.1976    
geneticcluster15 -0.160236   0.360150  -0.445   0.6588    
geneticcluster17 -0.160972   0.360150  -0.447   0.6574    
geneticcluster18 -0.005794   0.360150  -0.016   0.9872    
geneticcluster23 -0.138021   0.415865  -0.332   0.7417    
geneticcluster24 -0.488998   0.415865  -1.176   0.2468    
geneticcluster26 -0.307579   0.360150  -0.854   0.3983    
geneticcluster29 -0.067191   0.339553  -0.198   0.8442    
geneticcluster30 -0.341362   0.415865  -0.821   0.4167    
geneticcluster31  0.505684   0.339553   1.489   0.1445    
geneticcluster32  0.288553   0.415865   0.694   0.4919    
geneticcluster33  0.146502   0.322128   0.455   0.6518    
geneticcluster34 -0.242271   0.339553  -0.714   0.4798    
geneticcluster35 -0.102934   0.328770  -0.313   0.7559    
geneticcluster36  0.133531   0.415865   0.321   0.7499    
geneticcluster37 -0.423874   0.328770  -1.289   0.2049    
geneticcluster38 -0.387278   0.322128  -1.202   0.2365    
geneticcluster39  0.188641   0.360150   0.524   0.6034    
geneticcluster40 -0.206357   0.360150  -0.573   0.5699    
geneticcluster41 -0.234532   0.415865  -0.564   0.5760    
geneticcluster42 -0.366177   0.415865  -0.881   0.3840    
geneticcluster43 -0.317148   0.415865  -0.763   0.4503    
geneticcluster44  0.104091   0.328770   0.317   0.7532    
geneticcluster45 -0.024693   0.415865  -0.059   0.9530    
geneticcluster46  0.031292   0.328770   0.095   0.9247    
geneticcluster47  0.463011   0.415865   1.113   0.2724    
geneticcluster48  0.136576   0.415865   0.328   0.7444    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2941 on 39 degrees of freedom
Multiple R-squared:  0.6296,    Adjusted R-squared:  0.3256 
F-statistic: 2.071 on 32 and 39 DF,  p-value: 0.01551
Code 
set.seed(123)
model3.geneticcluster.sh<-lm(shannon_entropy~geneticcluster, data=metadata.pup)
summary(model3.geneticcluster.sh) #Adjusted R2 = 0.3214, p-value = 0.01667

Call:
lm(formula = shannon_entropy ~ geneticcluster, data = metadata.pup)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3184 -0.3026  0.0000  0.2512  1.4817 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       5.40745    0.74915   7.218 1.08e-08 ***
geneticcluster5  -0.90992    0.86505  -1.052   0.2993    
geneticcluster8  -0.76525    0.86505  -0.885   0.3818    
geneticcluster9  -0.32644    0.91752  -0.356   0.7239    
geneticcluster10  0.04159    1.05946   0.039   0.9689    
geneticcluster12  1.36980    0.91752   1.493   0.1435    
geneticcluster14 -0.34701    0.91752  -0.378   0.7073    
geneticcluster15 -0.05477    0.91752  -0.060   0.9527    
geneticcluster17 -0.33343    0.91752  -0.363   0.7183    
geneticcluster18  0.12337    0.91752   0.134   0.8937    
geneticcluster23 -1.22317    1.05946  -1.155   0.2553    
geneticcluster24 -0.79920    1.05946  -0.754   0.4552    
geneticcluster26 -0.55923    0.91752  -0.609   0.5457    
geneticcluster29  0.54477    0.86505   0.630   0.5325    
geneticcluster30 -0.24461    1.05946  -0.231   0.8186    
geneticcluster31  1.50061    0.86505   1.735   0.0907 .  
geneticcluster32  1.15868    1.05946   1.094   0.2808    
geneticcluster33  0.39073    0.82065   0.476   0.6366    
geneticcluster34 -0.81442    0.86505  -0.941   0.3523    
geneticcluster35  0.01257    0.83758   0.015   0.9881    
geneticcluster36  0.72209    1.05946   0.682   0.4995    
geneticcluster37 -0.84073    0.83758  -1.004   0.3217    
geneticcluster38 -0.21563    0.82065  -0.263   0.7941    
geneticcluster39  1.09131    0.91752   1.189   0.2415    
geneticcluster40 -0.22313    0.91752  -0.243   0.8091    
geneticcluster41 -1.61476    1.05946  -1.524   0.1355    
geneticcluster42 -0.86983    1.05946  -0.821   0.4166    
geneticcluster43 -0.27149    1.05946  -0.256   0.7991    
geneticcluster44  0.20563    0.83758   0.246   0.8074    
geneticcluster45  0.84185    1.05946   0.795   0.4317    
geneticcluster46  0.40958    0.83758   0.489   0.6276    
geneticcluster47  1.43795    1.05946   1.357   0.1825    
geneticcluster48  0.66626    1.05946   0.629   0.5331    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7492 on 39 degrees of freedom
Multiple R-squared:  0.6272,    Adjusted R-squared:  0.3214 
F-statistic: 2.051 on 32 and 39 DF,  p-value: 0.01667
Code 
set.seed(123)
model3.geneticcluster.pd<-lm(faith_pd~geneticcluster, data=metadata.pup)
summary(model3.geneticcluster.pd) #Adjusted R2 = 0.4175, p-value = 0.002517

Call:
lm(formula = faith_pd ~ geneticcluster, data = metadata.pup)

Residuals:
   Min     1Q Median     3Q    Max 
-4.431 -1.465  0.000  1.194  5.244 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       18.2616     2.6042   7.012 2.06e-08 ***
geneticcluster5   -4.5980     3.0071  -1.529   0.1343    
geneticcluster8   -1.3791     3.0071  -0.459   0.6491    
geneticcluster9   -1.8644     3.1895  -0.585   0.5622    
geneticcluster10  -1.2640     3.6830  -0.343   0.7333    
geneticcluster12   4.9958     3.1895   1.566   0.1254    
geneticcluster14  -4.0794     3.1895  -1.279   0.2085    
geneticcluster15  -1.0312     3.1895  -0.323   0.7482    
geneticcluster17  -0.7161     3.1895  -0.225   0.8235    
geneticcluster18   1.2232     3.1895   0.384   0.7034    
geneticcluster23   0.7106     3.6830   0.193   0.8480    
geneticcluster24  -2.4012     3.6830  -0.652   0.5182    
geneticcluster26  -1.5645     3.1895  -0.490   0.6265    
geneticcluster29  -0.1451     3.0071  -0.048   0.9618    
geneticcluster30  -2.9287     3.6830  -0.795   0.4313    
geneticcluster31   7.0062     3.0071   2.330   0.0251 *  
geneticcluster32   3.9851     3.6830   1.082   0.2859    
geneticcluster33   3.3048     2.8528   1.158   0.2537    
geneticcluster34  -0.5995     3.0071  -0.199   0.8430    
geneticcluster35   0.3302     2.9116   0.113   0.9103    
geneticcluster36   1.0913     3.6830   0.296   0.7686    
geneticcluster37  -1.5934     2.9116  -0.547   0.5873    
geneticcluster38  -2.6722     2.8528  -0.937   0.3547    
geneticcluster39   2.4435     3.1895   0.766   0.4482    
geneticcluster40  -0.8282     3.1895  -0.260   0.7965    
geneticcluster41   0.4892     3.6830   0.133   0.8950    
geneticcluster42  -3.0254     3.6830  -0.821   0.4164    
geneticcluster43  -2.0134     3.6830  -0.547   0.5877    
geneticcluster44   3.3906     2.9116   1.165   0.2513    
geneticcluster45   0.2377     3.6830   0.065   0.9489    
geneticcluster46   1.3387     2.9116   0.460   0.6482    
geneticcluster47   5.5730     3.6830   1.513   0.1383    
geneticcluster48   1.4064     3.6830   0.382   0.7046    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.604 on 39 degrees of freedom
Multiple R-squared:   0.68, Adjusted R-squared:  0.4175 
F-statistic:  2.59 on 32 and 39 DF,  p-value: 0.002517
Code 
#BETA DIVERSITY
set.seed(123)
bray_dist.pup = phyloseq::distance(psrt.pup, method="bray")
set.seed(123)
jacc_dist.pup = phyloseq::distance(psrt.pup, method="jaccard")
set.seed(123)
uu_dist.pup = phyloseq::distance(psrt.pup, method="uunifrac")
set.seed(123)
wu_dist.pup = phyloseq::distance(psrt.pup, method="wunifrac")
metadata.psrt.pup.gc <- as(sample_data(psrt.pup), "data.frame")
metadata.psrt.pup.gc$geneticcluster<-as.factor(metadata.psrt.pup.gc$geneticcluster)

set.seed(123)
pmod.pup.geneticcluster.bray<-adonis2(bray_dist.pup~geneticcluster, data=metadata.psrt.pup.gc, by="margin", perm=9999)
pmod.pup.geneticcluster.bray
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = bray_dist.pup ~ geneticcluster, data = metadata.psrt.pup.gc, permutations = 9999, by = "margin")
               Df SumOfSqs     R2     F Pr(>F)    
geneticcluster 32  13.5129 0.6821 2.615  1e-04 ***
Residual       39   6.2978 0.3179                 
Total          71  19.8107 1.0000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
pmod.pup.geneticcluster.jacc<-adonis2(jacc_dist.pup~geneticcluster, data=metadata.psrt.pup.gc, by="margin", perm=9999)
pmod.pup.geneticcluster.jacc
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = jacc_dist.pup ~ geneticcluster, data = metadata.psrt.pup.gc, permutations = 9999, by = "margin")
               Df SumOfSqs      R2      F Pr(>F)    
geneticcluster 32  15.6671 0.61137 1.9172  1e-04 ***
Residual       39   9.9593 0.38863                  
Total          71  25.6264 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
pmod.pup.geneticcluster.wu<-adonis2(wu_dist.pup~geneticcluster, data=metadata.psrt.pup.gc, by="margin", perm=9999)
pmod.pup.geneticcluster.wu
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = wu_dist.pup ~ geneticcluster, data = metadata.psrt.pup.gc, permutations = 9999, by = "margin")
               Df  SumOfSqs      R2      F Pr(>F)    
geneticcluster 32 0.0040358 0.65971 2.3628  1e-04 ***
Residual       39 0.0020817 0.34029                  
Total          71 0.0061175 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
set.seed(123)
pmod.pup.geneticcluster.uu<-adonis2(uu_dist.pup~geneticcluster, data=metadata.psrt.pup.gc, by="margin", perm=9999)
pmod.pup.geneticcluster.uu
Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 9999

adonis2(formula = uu_dist.pup ~ geneticcluster, data = metadata.psrt.pup.gc, permutations = 9999, by = "margin")
               Df SumOfSqs      R2      F Pr(>F)    
geneticcluster 32  11.2718 0.64608 2.2248  1e-04 ***
Residual       39   6.1747 0.35392                  
Total          71  17.4465 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Recreate Figure S1

Individual taxa bar plots displaying the bacterial families in the post-transfer gut microbiota in adult and juvenile bank voles arranged according to experimental group

Code 
physeq2.f.3 = filter_taxa(psrt.mom.time3, function(x) mean(x) > 0.1, TRUE)
physeq3.f.3 = transform_sample_counts(physeq2.f.3, function(x) x / sum(x) )
glom.f.3 <- tax_glom(physeq3.f.3, taxrank = 'Family')
data_glom.3 <- psmelt(glom.f.3) 
data_glom.3$Family <- as.character(data_glom.3$Family) 
data_glom.3$Family <- factor(data_glom.3$Family, 
                             levels = c("NA","Campylobacteraceae",
                                        "Ruminococcaceae", "Christensenellaceae","Desulfovibrionaceae",
                                        "Oscillospiraceae","Lachnospiraceae",
                                        "Lactobacillaceae","Muribaculaceae"))

data_glom.3$treatment[data_glom.3$treatment=="urburb"]<-"Urban-Urban"
data_glom.3$treatment[data_glom.3$treatment=="urbcont"]<-"Urban-Rural"
data_glom.3$treatment[data_glom.3$treatment=="conturb"]<-"Rural-Urban"
data_glom.3$treatment[data_glom.3$treatment=="contcont"]<-"Rural-Rural"
plot.fig.S3<-ggplot(data_glom.3, 
       aes(x=Sample, y=Abundance, fill=Family)) + 
  geom_bar(stat="identity", position="fill",color="black") + 
  facet_grid(~treatment, scale="free")+
  scale_fill_manual(values = c("#9E0142", "#D53E4F", "#F46D43", "#FDAE61", "#FEE08B","#ABDDA4","blue","purple")) +
  theme_bw()+
  xlab("Individual adult bank vole")+
  ylab("Relative Abundance")+
  theme(axis.text.x=element_blank())

physeq2.f.pup = filter_taxa(psrt.pup, function(x) mean(x) > 0.1, TRUE)
physeq3.f.pup = transform_sample_counts(physeq2.f.pup, function(x) x / sum(x) )
glom.f.pup <- tax_glom(physeq3.f.pup, taxrank = 'Family')
data_glom.pup <- psmelt(glom.f.pup) 
data_glom.pup$Family <- as.character(data_glom.pup$Family) 
data_glom.pup$Family <- factor(data_glom.pup$Family, 
                             levels = c("NA","Campylobacteraceae",
                                        "Ruminococcaceae", "Christensenellaceae","Desulfovibrionaceae",
                                        "Oscillospiraceae","Lachnospiraceae",
                                        "Lactobacillaceae","Muribaculaceae"))

data_glom.pup$treatment[data_glom.pup$treatment=="urburb"]<-"Urban-Urban"
data_glom.pup$treatment[data_glom.pup$treatment=="urbcont"]<-"Urban-Rural"
data_glom.pup$treatment[data_glom.pup$treatment=="conturb"]<-"Rural-Urban"
data_glom.pup$treatment[data_glom.pup$treatment=="contcont"]<-"Rural-Rural"
plot.fig.S3.pup<-ggplot(data_glom.pup, 
                    aes(x=Sample, y=Abundance, fill=Family)) + 
  geom_bar(stat="identity", position="fill",color="black") + 
  facet_grid(~treatment, scale="free")+
  scale_fill_manual(values = c("#9E0142", "#D53E4F", "#F46D43", "#FDAE61", "#FEE08B","#ABDDA4","blue","purple")) +
  theme_bw()+
  xlab("Individual pup")+
  ylab("Relative Abundance")+
  theme(axis.text.x=element_blank())

plot.fig.S3.mom.pup<-ggarrange(plot.fig.S3, plot.fig.S3.pup,
                       common.legend = TRUE,
                       legend = "right",
                       labels = c("A", "B"),
                       ncol = 1, nrow = 2)

plot.fig.S3.mom.pup # 8 x 6

Recreate Figure S2

PCoA ordinations visualising the compositions of the pre-transfer gut microbiota of adult bank voles

Code 
p2.ord.bray.mom.1 <- ordinate(psrt.mom.time1, "PCoA", distance=bray_dist.mom.time1) 
bray.mom.1 = plot_ordination(psrt.mom.time1, p2.ord.bray.mom.1, type="samples", color="sitetr1full",  axes=c(1,2)) +
  stat_ellipse(aes(color=sitetr1), alpha=.5,linewidth =2)+
  scale_color_manual(values=c("#693aa3","green", "#f00078"), name=c(""),
                     labels=c("Rural (Jyväskylä)","Rural (Konnevesi)", "Urban"))+
  xlab("PCoA Axis.1  [16%]")+
  ylab("PCoA Axis.2  [9.6%]")+
  ggtitle("Bray-Curtis: Pre-transfer") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

p2.ord.jacc.mom.1 <- ordinate(psrt.mom.time1, "PCoA", distance=jacc_dist.mom.time1) 
jacc.mom.1 = plot_ordination(psrt.mom.time1, p2.ord.jacc.mom.1, type="samples", color="sitetr1full",  axes=c(1,2)) +
  stat_ellipse(aes(color=sitetr1), alpha=.5,linewidth =2)+
  scale_color_manual(values=c("#693aa3","green", "#f00078"), name=c(""),
                     labels=c("Rural (Jyväskylä)","Rural (Konnevesi)", "Urban"))+
  xlab("PCoA Axis.1  [11%]")+
  ylab("PCoA Axis.2  [7.2%]")+
  ggtitle("Jaccard: Pre-transfer") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

p2.ord.wunifrac.mom.1 <- ordinate(psrt.mom.time1, "PCoA", distance=wunifrac_dist.mom.time1) 
wunifrac.mom.1 = plot_ordination(psrt.mom.time1, p2.ord.wunifrac.mom.1, type="samples", color="sitetr1full",  axes=c(1,2)) +
  stat_ellipse(aes(color=sitetr1), alpha=.5,linewidth =2)+
  scale_color_manual(values=c("#693aa3","green", "#f00078"), name=c(""),
                     labels=c("Rural (Jyväskylä)","Rural (Konnevesi)", "Urban"))+
  xlab("PCoA Axis.1  [45.2%]")+
  ylab("PCoA Axis.2  [12%]")+
  ggtitle("wUnifrac: Pre-transfer") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

p2.ord.uunifrac.mom.1 <- ordinate(psrt.mom.time1, "PCoA", distance=uunifrac_dist.mom.time1) 
uunifrac.mom.1 = plot_ordination(psrt.mom.time1, p2.ord.uunifrac.mom.1, type="samples", color="sitetr1full",  axes=c(1,2)) +
  stat_ellipse(aes(color=sitetr1), alpha=.5,linewidth =2)+
  scale_color_manual(values=c("#693aa3","green", "#f00078"), name=c(""),
                     labels=c("Rural (Jyväskylä)","Rural (Konnevesi)", "Urban"))+
  xlab("PCoA Axis.1  [10.2%]")+
  ylab("PCoA Axis.2  [9.8%]")+
  ggtitle("unwUnifrac: Pre-transfer") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

plot.fig.t1<-ggarrange(bray.mom.1, jacc.mom.1,wunifrac.mom.1,uunifrac.mom.1,
                       common.legend = TRUE,
                       legend = "right",
                       labels = c("A", "B","C","D"),
                       ncol = 2, nrow = 2)

plot.fig.t1 # 7 x 6

Recreate Figure S3

Change in alpha diversity of the gut microbiota of bank voles throughout the experiment (ASV Richness)

Code 
plot.asv.t1.to.t3<-ggplot(data=diff.13, aes(x=treatment, y=asv_1_to_3, fill=treatment))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "ASV richness", x= "",y = "Change in alpha diversity")+
  scale_x_discrete(labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  scale_fill_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"),name=c("RT exp. group"),
                    labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  theme_minimal()

plot.asv.t3<-ggplot(data=metadata.mom.time3, aes(x=treatment, y=observed_features, fill=treatment))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "Post-transfer ASV richness", x= "",y = "Total alpha diversity")+
  scale_x_discrete(labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  scale_fill_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"),name=c("RT exp. group"),
                    labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  theme_minimal()

plot.alpha.change<-ggarrange(plot.asv.t1.to.t3+ rremove("x.text"), plot.asv.t3 + rremove("x.text"),
                             common.legend = TRUE,
                             legend = "right",
                             labels = c("A", "B"),
                             ncol = 2, nrow = 1)

plot.alpha.change

Recreate Figure S4

Magnitude of change in beta diversity of the gut microbiota of bank voles in terms of past habitat

Code 
plot.bray.t1.to.t3<-ggplot(data=dist.13, aes(x=sitetr1, y=bray_1_to_3, fill=sitetr1))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "Bray-Curtis metric", x= "",y = "Change in Beta Diversity")+
  scale_x_discrete(labels=c("Rural","Urban"))+
  scale_fill_manual(values=c("#693aa3","#f00078"),name=c(""),
                    labels=c("Rural","Urban"))+
  geom_signif(y_position = c(0.95), xmin = c(1), xmax = c(2),annotation = c("ns"), tip_length = 0)+
  theme_minimal()

plot.jacc.t1.to.t3<-ggplot(data=dist.13, aes(x=sitetr1, y=jacc_1_to_3, fill=sitetr1))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "Jaccard metric", x= "",y = "Change in Beta Diversity")+
  scale_x_discrete(labels=c("Rural","Urban"))+
  scale_fill_manual(values=c("#693aa3","#f00078"),name=c(""),
                    labels=c("Rural","Urban"))+
  geom_signif(y_position = c(0.875), xmin = c(1), xmax = c(2),annotation = c("ns"), tip_length = 0)+
  theme_minimal()

plot.wunifrac.t1.to.t3<-ggplot(data=dist.13, aes(x=sitetr1, y=wunifrac_1_to_3, fill=sitetr1))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "wUnifrac metric", x= "",y = "Change in Beta Diversity")+
  scale_x_discrete(labels=c("Rural","Urban"))+
  scale_fill_manual(values=c("#693aa3","#f00078"),name=c(""),
                    labels=c("Rural","Urban"))+
  geom_signif(y_position = c(0.45), xmin = c(1), xmax = c(2),annotation = c("*"), tip_length = 0)+
  theme_minimal()

plot.unwunifrac.t1.to.t3<-ggplot(data=dist.13, aes(x=sitetr1, y=unwunifrac_1_to_3, fill=sitetr1))+
  geom_boxplot() +
  theme(legend.position="right")+
  labs(title= "unwUnifrac metric", x= "",y = "Change in Beta Diversity")+
  scale_x_discrete(labels=c("Rural","Urban"))+
  scale_fill_manual(values=c("#693aa3","#f00078"),name=c(""),
                    labels=c("Rural","Urban"))+
  geom_signif(y_position = c(0.575), xmin = c(1), xmax = c(2),annotation = c("ns"), tip_length = 0)+
  theme_minimal()

plot.fig.turnover<-ggarrange(plot.bray.t1.to.t3+ rremove("x.text"), plot.jacc.t1.to.t3 + rremove("x.text"),plot.wunifrac.t1.to.t3+ rremove("x.text"),plot.unwunifrac.t1.to.t3+ rremove("x.text"),
                             common.legend = TRUE,
                             legend = "right",
                             labels = c("A", "B","C","D"),
                             ncol = 2, nrow = 2)

plot.fig.turnover

Recreate Figure S5

Additional CAP ordinations visualising the compositions of adult and juvenile post-transfer gut microbiota in terms of past and present habitat

Code 
#######
#ADULTS
#######
set.seed(123)
cap.mom.3.jacc <- ordinate(physeq = psrt.mom.time3, method = "CAP",distance = jacc_dist.mom.time3,formula = ~ sitetr1+sitetr2)
plot.mom.3.jacc<-plot_ordination(psrt.mom.time3, cap.mom.3.jacc,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Adults: Jaccard metric")+xlab("CAP Axis.1  [6.4%]")+ylab("CAP Axis.2  [3.7%]")+theme_minimal()

set.seed(123)
cap.mom.3.wunifrac <- ordinate(physeq = psrt.mom.time3, method = "CAP",distance = wunifrac_dist.mom.time3,formula = ~ sitetr1+sitetr2)
plot.mom.3.wunifrac<-plot_ordination(psrt.mom.time3, cap.mom.3.wunifrac,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Adults: wUnifrac metric")+xlab("CAP Axis.1  [7.9%]")+ylab("CAP Axis.2  [4.5%]")+theme_minimal()

set.seed(123)
cap.mom.3.unwunifrac <- ordinate(physeq = psrt.mom.time3, method = "CAP",distance = uunifrac_dist.mom.time3,formula = ~ sitetr1+sitetr2)
plot.mom.3.uunifrac<-plot_ordination(psrt.mom.time3, cap.mom.3.unwunifrac,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Adults: unwUnifrac metric")+xlab("CAP Axis.1  [4.8%]")+ylab("CAP Axis.2  [4.4%]")+theme_minimal()

plot.cap.adults.si<-ggarrange(plot.mom.3.jacc, plot.mom.3.wunifrac,plot.mom.3.uunifrac,
                              common.legend = TRUE,
                              legend = "right",
                              labels = c("A","B", "C"),
                              ncol = 3, nrow = 1, align="v"+ theme_void())

########
#OFFSPRING
########
set.seed(123)
jacc_dist.pup.all = phyloseq::distance(psrt.pup, method="jaccard")
set.seed(123)
wunifrac_dist.pup.all = phyloseq::distance(psrt.pup, method="wunifrac")
set.seed(123)
uunifrac_dist.pup.all = phyloseq::distance(psrt.pup, method="uunifrac")
metadata.psrt.pup.all <- as(sample_data(psrt.pup), "data.frame")
set.seed(123)
cap.pup.all.jacc <- ordinate(physeq = psrt.pup, method = "CAP",distance = jacc_dist.pup.all,formula = ~ sitetr1+sitetr2)
set.seed(123)
cap.pup.all.wunifrac <- ordinate(physeq = psrt.pup, method = "CAP",distance = wunifrac_dist.pup.all,formula = ~ sitetr1+sitetr2)
set.seed(123)
cap.pup.all.unwunifrac <- ordinate(physeq = psrt.pup, method = "CAP",distance = uunifrac_dist.pup.all,formula = ~ sitetr1+sitetr2)

plot.pup.cap.jacc<-plot_ordination(psrt.pup, cap.pup.all.jacc,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Offspring: Jaccard metric")+xlab("CAP Axis.1  [3%]")+ylab("CAP Axis.2  [1.8%]")+theme_minimal()

plot.pup.cap.wunifrac<-plot_ordination(psrt.pup, cap.pup.all.wunifrac,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Offspring: wUnifrac metric")+xlab("CAP Axis.1  [2.9%]")+ylab("CAP Axis.2  [1.9%]")+theme_minimal()

plot.pup.cap.unwunifrac<-plot_ordination(psrt.pup, cap.pup.all.unwunifrac,type = "samples",color="treatment")+geom_point(size = 3)+
  scale_color_manual(values=c("#693aa3","#7777ff","#ffc1e0","#f00078"), name=(""),
                     labels=c("Rural-Rural","Rural-Urban","Urban-Rural","Urban-Urban"))+
  ggtitle("Offspring: unwUnifrac metric")+xlab("CAP Axis.1  [3%]")+ylab("CAP Axis.2  [2.5%]")+theme_minimal()

plot.cap.offspring.si<-ggarrange(plot.pup.cap.jacc, plot.pup.cap.wunifrac,plot.pup.cap.unwunifrac,
                              common.legend = TRUE,
                              legend = "right",
                              labels = c("D","E", "F"),
                              ncol = 3, nrow = 1, align="v"+ theme_void())

plot.cap.all.si<-ggarrange(plot.cap.adults.si, plot.cap.offspring.si,
                              common.legend = TRUE,
                              legend = "right",
                              ncol = 1, nrow = 2, align="v"+ theme_void())
plot.cap.all.si # 10 x 12

Recreate Figure S6

PCoA ordinations visualising the compositions of adult and juvenile post-transfer gut microbiota in terms of genetic cluster

Code 
set.seed(123)
bray_dist.pup = phyloseq::distance(psrt.pup, method="bray")
set.seed(123)
jacc_dist.pup = phyloseq::distance(psrt.pup, method="jaccard")
set.seed(123)
uu_dist.pup = phyloseq::distance(psrt.pup, method="uunifrac")
set.seed(123)
wu_dist.pup = phyloseq::distance(psrt.pup, method="wunifrac")

#MATCHED mothers + offspring
set.seed(123)
bray_dist.matched = phyloseq::distance(psrt.matched, method="bray")
set.seed(123)
jacc_dist.matched = phyloseq::distance(psrt.matched, method="jaccard")
set.seed(123)
uunifrac_dist.matched = phyloseq::distance(psrt.matched, method="uunifrac")
set.seed(123)
wunifrac_dist.matched = phyloseq::distance(psrt.matched, method="wunifrac")

p2.ord.bray.matched <- ordinate(psrt.matched, "PCoA", distance=bray_dist.matched) 
bray.matched = plot_ordination(psrt.matched, p2.ord.bray.matched, type="samples", color="geneticcluster_2",  shape="type",axes=c(1,2)) +
   geom_point(alpha=2, size=3)+
  labs(color='Genetic cluster',shape="Type")+ 
  xlab("PCoA Axis.1  [15.5%]")+
  ylab("PCoA Axis.2  [10.7%]")+
  ggtitle("Matched mothers and offspring (Bray-Curtis metric)") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

p2.ord.bray.pup <- ordinate(psrt.pup, "PCoA", distance=bray_dist.pup) 
plot.pcoa.pup.bray = plot_ordination(psrt.pup, p2.ord.bray.pup, type="samples", color="geneticcluster_2",axes=c(1,2)) +
  geom_point(alpha=2, size=3)+
  labs(color='Genetic cluster')+ 
  xlab("PCoA Axis.1  [11.3%]")+
  ylab("PCoA Axis.2  [9.5%]")+
  ggtitle("All offspring (Bray-Curtis metric)") +
  scale_y_continuous(position = "left")+
  theme_minimal() 

plot.geneticcluster<-ggarrange(bray.matched, plot.pcoa.pup.bray,
                              common.legend = FALSE,
                              legend = "right",
                              labels = c("A","B"),
                              ncol = 2, nrow = 1, align="v"+ theme_void())
plot.geneticcluster # 10 x 12